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THE  APPLICATION  OF  SECULAR  PERTURBATION  THEORY  TO 
EXPLAIN  WARPING  IN  THE  CIRCUMSTELLAR  DISK  OF  BETA  PICTORIS 

By 

Steven  J.  Novotny 
December  2003 

Chair:  Stanley  F.  Dermott 
Major  Department:  Astronomy 

This  research  is  a  numerical  investigation  into  the  dynamical  influences  of  planets  on 
the  dust  disks  surrounding  young  main  sequence  stars.  Motivating  this  research  effort  are 
the  observations  of  the  Beta  Pictoris  circumstellar  disk  made  by  the  UF  Department  of 
Astronomy’s  mid-IR  team.  These  IR  images  show  with  unprecedented  clarity  the  features 
and  asymmetries  of  the  inner,  <  100  AU,  portion  of  the  Beta  Pictoris  disk;  the  most 
interesting  of  which  is  the  dramatic  warping  of  the  disk’s  mid-plane.  Analyses  of  prior 
observations  have  suggested  that  the  features  are  attributed  to  the  presence  of  a  planet 
or  a  planetary  system.  Past  dynamic  analysis  has  focused  primarily  on  the  presence  of  a 
single  planet  and  the  resulting  perturbations  on  the  dust  disk  through  hydrodynamic  or 
N-body  analysis.  This  research  will  show  that  the  type  of  features  observed  in  these  images, 
specifically  the  warping  of  the  disk,  can  also  be  explained  with  a  system  of  two  (or  more) 
planets  and  secular  perturbation  theory  while  using  more  plausible  assumptions  than  did 
other  models. 
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CHAPTER  1 
INTRODUCTION 


1.1  Motivation 

As  of  18  December,  2002,  the  discoveries  of  102  exosolar  planets  have  been  claimed 
(Schneider,  2002).  These  planets  exist  in  88  planetary  systems,  11  of  which  are  multiple 
planet  systems.  All  have  been  detected  through  careful  observations  of  dynamical  effects 
with  the  vast  majority  detected  through  stellar  radial  velocity  measurements  (Perryman, 
2002).  Radial  velocity  measurements  are  made  by  observing  the  doppler  shift  in  the  radiated 
light  from  a  star.  The  shifts  in  wavelength  result  from  the  gravitational  perturbation  of 
a  planet  on  the  star  around  which  it  orbits.  Using  indirect  methods,  such  as  observing 
the  effect  of  a  planet  on  its  central  star  instead  of  the  planet  itself,  is  necessary  as  current 
technology  has  not  yet  allowed  direct  observation  of  an  exosolar  planet.  Current  direct 
detection  methods  with  the  necessary  resolution  and  sensitivity  are  not  yet  in  operation. 

Several  projects  are  underway  to  develop  the  capability  to  make  direct  detections  of  ex¬ 
osolar  planets.  The  most  significant  of  these  projects  include  NASA’s  Space  Interferometry 
Mission  (SIM),  the  Terrestrial  Planet  Finder  (TPF),  and  ESAs  Infrared  Space  Interferome¬ 
try  Mission  (IRSI).  All  of  these  systems  will  use  interferometry  to  achieve  the  high  angular 
resolution  needed  for  planetary  detection.  Interferometry  uses  the  interference  of  light  as¬ 
sociated  with  the  observations  from  two  or  more  different  detectors.  The  goal  of  designing 
an  interferometer  is  to  have  the  greatest  possible  baseline  (the  distance  between  the  compo¬ 
nent  detectors)  as  the  angular  resolution  achieved  with  the  instrument  is  a  function  of  the 
baseline  distance.  SIM  is  a  space  based  optical  interferometer  with  a  ten  meter  baseline  and 
is  scheduled  for  launch  in  2009.  The  final  design  for  TPF  will  be  decided  in  2006  and  the 
launch  date  is  somewhere  between  2012  and  2015.  The  IRSI  is  the  most  ambitious  project 
and  currently  plans  to  utilize  six  different  space  telescopes  providing  an  effective  baseline 
of  60  meters  and  is  expected  to  launch  around  2014.  Though  these  projects  will  provide 


1 


2 


great  capability,  they  are  all  still  far  from  being  operational.  At  this  point,  the  observation 
of  individual  planets  is  still  out  of  reach  with  current  instruments,  and  thus  the  indirect 
methods  such  as  radial  velocity  measurements  and  observing  brightness  variations  are  still 
required. 

Recent  observations  of  circumstellar  disks  have  provided  another  possibility  for  plane¬ 
tary  detection.  One  example  is  the  disk  surrounding  the  star  HR4796A  which  was  resolved 
in  the  mid-IR  (Telesco  et  al.,  2000).  The  nearly  edge-on  view  of  this  disk  showed  a  double 
lobed  feature  that  exhibited  a  possible  overall  asymmetry  in  brightness  between  the  two 
lobes  (Wyatt  et  al.,  1999).  This  asymmetry  was  interpreted  as  the  disk  being  offset  from 
the  central  star.  As  will  be  discussed  later,  the  measurement  of  this  asymmetry  was  used 
to  make  assertions  about  the  possibility  of  the  disk  being  influenced  by  the  presence  of  a 
planet.  Another  example  is  Beta  Pictoris  -  the  most  well  studied  of  all  circumstellar  disks. 
Beta  Pictoris  is  an  approximately  20  Myr  old  main  sequence  star  with  a  disk  that  is  large 
(~  1000  AU),  bright,  and  viewed  nearly  edge-on.  Even  more  interesting  than  the  size  and 
brightness  of  the  Beta  Pictoris  disk  are  the  large  and  small-scale  asymmetries  and  features 
observed  in  the  morphology  of  the  disk  -  the  most  interesting  being  the  warping  of  the 
inner  disk. 

The  mid-plane  warping  of  the  disk  is  interesting  as  it  is  arguably  the  distinct  signature 
of  gravitational  perturbations.  Several  models  have  been  developed  to  explain  the  source 
of  these  perturbations.  The  most  in-depth  work  has  been  based  on  the  assumption  of  a 
close  encounter  with  another  star  (Kalas  et  al.,  2000),  or  the  presence  of  a  single  planet 
on  an  orbit  inclined  to  the  plane  of  the  disk  (Augereau  et  al.,  2001).  These  models  do  in 
fact  account  for  warping  in  the  model  disks,  but  there  are  several  unresolved  issues.  First, 
there  is  no  observed  star  that  meets  the  necessary  criteria  of  proximity  and  proper  motion 
to  match  the  fly-by  star.  Second,  there  is  no  clear  explanation  as  to  why  a  single  planet 
would  be  in  an  orbit  with  a  significant  inclination  to  the  disk  from  which  it  formed.  Current 
theories  of  planet  formation  assert  that  planets  form  from  the  disk  material  and  therefore, 
without  outside  perturbations,  should  be  nearly  coplanar  with  the  disk. 

An  alternative  approach  in  investigating  the  source  of  the  warping  is  based  on  a  mech¬ 
anism  that  explains  the  observed  warping  in  the  zodiacal  cloud  of  our  own  solar  system 
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-  secular  perturbations  from  the  planets.  Secular  perturbations  are  long  term  effects  that 
are  the  result  of  the  combined  gravitational  perturbations  of  planets.  Secular  perturbations 
manifest  themselves  as  a  contribution  to  the  osculating  orbital  elements  of  the  individual 
dust  particles  that  comprise  the  disk.  More  simply  stated,  this  approach  is  based  on  the 
premise  that  the  cumulative  effects  of  gravitational  perturbations,  on  long  timescales,  will 
play  a  significant  role  in  determining  the  observed  characteristics  of  a  dust  particle’s  orbital 
trajectory,  and  hence  the  morphology  of  the  disk.  In  this  research  a  second  order  expan¬ 
sion,  in  inclination  and  eccentricity,  of  the  distmbing  function  is  used  to  calculate  these 
forced  components.  The  theoretical  framework  for  determining  these  values  is  provided 
and  discussed  in  chapter  3. 

To  support  the  theory  of  secular  perturbations  being  the  dominant  influence  on  the 
morphology  of  a  disk,  such  as  that  of  Beta  Pictoris,  it  necessary  to  show  that  the  conditions 
within  the  disk  are  conducive  to  this  process.  For  this  research,  this  is  done  by  determining 
the  timescales  associated  with  each  of  the  physical  processes  expected  to  be  occurring  in 
a  disk  with  the  attributes  observed  for  Beta  Pictoris.  Some  of  the  physical  processes  con¬ 
sidered  include  the  Poynting-Robertson  drag,  collisional  destruction,  and  gas  drag.  These 
processes  axe  dependent  on  characteristics  of  the  star/disk  system  such  as  the  star’s  lumi¬ 
nosity,  the  particle  density  within  the  disk,  and  the  gas  content  respectively.  These  effects 
and  the  associated  analyses  are  explained  further  in  chapter  2. 

The  following  section  provides  brief  examples  of  observed  asymmetries  in  circumstellar 
disks  and  discusses  how  secular  perturbation  theory  can  be  used  to  explain  the  observed 
features. 


1.2  Observed  Asymmetries  in  Circumstellar  Disks 
Indirect  methods  for  determining  the  presence  of  a  planetary  system  are  still,  by  far,  the 
most  commonly  used  in  the  search  for  exosolar  planets.  However,  when  we  are  able  to  resolve 
the  disk  and  characterize  the  details  of  the  morphology,  we  have  the  rare  opportunity  to 
probe  the  dynamics  of  the  material  within  the  disk  and  the  dynamical  influences  of  a  planet 
or  planetary  system.  Two  such  examples  are  the  disks  smrounding  the  stars  HR4796A  and 
Beta  Pictoris.  These  systems  are  well  suited  for  the  application  of  secular  perturbation 
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theory  as  the  disks  surrounding  theses  stars  have  arguably  existed  long  enough  for  secular 
effects  to  have  influenced  the  entire  disk. 

When  looking  for  dynamical  influences  on  secular  timescales,  there  are  two  key  observ¬ 
able  features:  pericenter  glow  and  mid-plane  warping.  As  will  be  discussed  next,  pericenter 
glow  is  the  result  of  a  disk  being  offset  from  the  star  due  to  the  eccentricity  imposed  upon 
the  orbits  of  the  disk’s  dust  particles  by  a  perturbing  planet.  Similarly,  warping  will  be 
shown  to  be  the  result  of  an  imposed  inclination. 

Essential  for  understanding  secular  perturbation  theory  is  the  idea  that  the  observable, 
or  osculating,  orbital  elements  of  a  dust  particle  within  a  disk  are  a  combination  of  forced 
and  proper  elements.  The  forced  elements  are  associated  with  the  gravitational  perturba¬ 
tions  from  the  planets  and  the  proper  elements  are  associated  with  the  trajectory  of  the 
particles  at  the  time  they  were  created.  This  will  be  described  in  detail  in  chapter  3. 

1.2.1  Pericenter  Glow 

The  phenomenon  of  pericenter  glow  was  investigated  by  Wyatt  et  al.  (1999)  in  an  effort 
to  explain  the  brightness  asymmetry  in  HR4796.  This  work  was  motivated  by  Telesco  et  al. 
(1999)  who  identifled  the  5.9  ±  3.2%  brightness  asymmetry  in  the  double-lobed  structure  of 
the  disk.  This  double-lobed  structure  was  considered  to  be  the  result  of  viewing  the  disk 
nearly  edge  on  -  tilted  approximately  13  degrees  to  the  line  of  sight. 

Wyatt  et  al.  proposed  that  the  secular  perturbations  of  a  single  massive  planet  imposed 
a  forced  eccentricity  upon  the  individual  particle  orbits,  such  that  the  overall  disk’s  center 
of  symmetry  was  offset  from  the  central  star.  This  offset  resulted  in  one  portion  of  the 
disk  being  closer  to  the  star,  and  thus  hotter  and  brighter.  This  difference  resulted  in  the 
observed  brightness  asymmetry.  Detailed  modeling  showed  that  a  forced  eccentricity  as 
small  as  0.02  could  well  account  for  the  ~  5%  lobe  brightness  asymmetry. 

Another  important  conclusion  from  the  work  of  Wyatt  et  al.  is  that  the  observed  dust 
particles  are  actually  a  tracer  of  the  larger  bodies.  By  modeling  the  dust  particles  as  Mie 
spheres,  it  was  determined  that  the  particles  responsible  for  the  mid-IR  flux  typically  have 
diameters  of  2  —  3  pm.  However,  these  particles  would  be  quickly  removed  from  the  system 
due  to  radiation  pressure  and  therefore  cannot  be  primordial.  This  implies  a  source  for  the 


5 


smaller  particles,  arguably  large  bodies  whose  dynamics  are  dominated  by  secular  pertur¬ 
bations  over  radiation  based  forces. 

1.2.2  Mid-plane  Warping 

The  best  example  of  a  directly  observed  warp  in  a  circmnstellar  disk  is  in  that  of  Beta 
Pictoris.  Since  its  initial  IR  imaging  by  Lagage  and  Pantin  (1994),  warping  of  the  inner  disk 
(<  100  AU)  was  apparent.  One  of  the  best  characterizations  of  the  mid-plane  warping  was 
accomplished  by  Heap  et  al.  (2000)  who  imaged  the  disk  with  the  Space  Telescope  Infrared 
Spectrograph  (STIS)  on  the  Hubble  Space  Telescope.  Their  observations  indicated  warping 
that  extended  out  to  ~80  AU  and  reached  a  maximum  amplitude  of  ~1.5  AU.  In  addition 
to  clearly  showing  the  warping,  they  quantified  the  warp  by  illustrating  how  it  could  be 
well  modeled  as  the  superposition  of  two  concentric  disk-like  components,  one  inclined  to 
the  other  by  rsj  3  degrees.  As  will  be  discussed  later,  several  dynamical  explanations  were 
proposed  to  explain  the  two  component  model. 

More  recent  imagiTig  of  Beta  Pictoris  (Weinberger  et  al.  2003,  Wahhaj  et  al.  2003) 
indicated  additional  warping  of  the  inner  region  of  the  disk  (<  20  AU)  which  was  not 
observed  by  Heap  et  al.  Furthermore,  the  inner  disk  warping  was  opposite  in  direction 
to  the  large  scale  warping  observed  in  the  STIS  images.  As  will  be  discussed  later,  this 
additional  warping  further  complicates  the  models  based  on  a  single  planet. 

Within  this  work  it  will  be  argued  that  mid-plane  warping  is  the  result  of  a  mechanism 
similar  to  that  which  causes  pericenter  glow.  Whereas  an  imposed  forced  eccentricity  causes 
the  disk  to  be  offset  from  the  central  star,  an  imposed  forced  inclination  causes  warping. 
This  phenomenon  has  been  observed  and  studied  for  the  solar  system’s  zodiacal  cloud  (Der- 
mott  et  al.,  2001).  The  general  idea  behind  warping  is  that  a  system  consisting  of  multiple 
planets  will  impose  a  forced  inclination  upon  the  orbits  of  dust  particles  in  a  way  that  is 
dependent  on  semi-major  axis.  Furthermore,  the  forced  inclination  defines  the  planes  about 
which  the  perturbed  particle  orbits  process.  Since  the  forced  inclination  defines  the  dust 
particles’  mean  orbital  planes  and  since  the  forced  inclination  varies  with  semi-major  axis, 
warping  occurs.  Determination  of  the  forced  orbital  elements  arising  from  planetary  per¬ 
turbations  is  accomplished  through  secular  perturbation  theory.  The  application  of  secular 
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perturbation  theory  to  the  mid-plane  warping  of  circumstellar  disks  is  discussed  in  detail 
in  the  following  chapters. 


1.3  Dissertation  Overview 

The  remainder  of  this  dissertation  will  discuss  the  physics  associated  with  the  study  of 
the  circumstellar  environment  and  application  of  such  models  in  interpreting  the  observed 
characteristics  of  a  complex  system  such  as  the  disk  of  Beta  Pictoris.  Chapter  2  reviews  the 
current  theories  of  circumstellar  disk  formation  and  the  subsequent  planet  formation.  Also 
provided  within  this  chapter,  is  a  review  of  the  physical  processes  that  affect  the  distribution 
of  dust  within  the  disk.  Chapter  3  will  discuss  the  theoretical  framework  for  determining 
the  long-term  effects  of  gravity  on  dust  particle  trajectories.  This  is  done  through  a  review 
of  secular  perturbation  theory  using  a  second  order  expansion  of  the  disturbing  function  in 
eccentricity  and  inclination.  The  resulting  solution  in  terms  of  the  classical  orbital  elements 
is  presented.  Chapter  4  explains  the  implications  of  these  results  with  regard  to  the  effects 
on  disk  morphology.  Special  emphasis  is  placed  on  the  role  that  a  forced  inclination  that 
varies  with  semi-major  axis  plays  in  determining  the  warping  of  a  disk’s  mid-plane.  Chapter 
5  reviews  the  observational  history  of  Beta  Pictoris  starting  with  the  initial  detection  of  the 
disk  from  the  initial  IRAS  measurements  of  the  spectral  energy  distribution  to  the  most 
recent  high  resolution  images  of  the  disk  made  with  the  Hubble  Space  Telescope.  Also 
discussed  in  this  chapter  are  the  theoretical  models  used  to  explain  the  unique  features  and 
asymmetries  indicated  in  these  observations. 

The  final  three  chapters,  chapters  7,  8,  and  9,  present  the  application  of  the  theoretical 
tools  developed  in  previous  chapters  to  address  the  specifics  of  observed  features  -  not  only 
in  the  prior  observations,  but  in  the  most  recent  high  resolution  mid-IR  images  of  Beta 
Pictoris.  Also  addressed  are  limitations  to  the  models  developed  to  simulate  the  observed 
dust  disk.  These  issues  include  consideration  of  a  massive  disk  and  the  possible  effects 
associated  with  a  non-spherical,  or  oblate,  central  star. 


CHAPTER  2 

THE  CIRCUMSTELLAR  ENVIRONMENT 


2.1  Introduction 

The  circumstellar  environment  changes  greatly  over  the  lifetime  of  the  star.  In  par¬ 
ticular,  the  cloud  of  material  surrounding  stars,  mostly  gas  and  dust,  undergoes  great 
transformations  in  both  composition  and  spatial  distribution.  Observational  evidence  sug¬ 
gests  that  as  a  normal  result  of  the  star  formation  process  that  all  stars  may  possibly  be 
born  with  a  distribution  of  dust  and  gas  in  the  form  of  a  disk  (Beckwith,  1999).  However, 
over  time  the  material  is  subjected  to  many  influences  that  may  result  in  the  partial  or 
complete  destruction  of  a  well  structured,  easily  observable  disk.  One  example  of  the  many 
changes  that  occur  is  in  the  composition  of  the  disk.  When  a  protostar  reaches  the  main 
sequence  it  is  estimated  that  all  but  one  percent  of  the  mass  of  the  disk  is  in  hydrogen 
and  helium  gas.  The  remaining  one  percent  is  in  the  form  of  dust  that  is  closely  coupled 
with  the  gas.  Over  time,  the  gas  is  removed  and  the  dust  continues  to  evolve.  The  dust 
particles,  depending  on  size,  are  generally  destroyed;  blown  away  by  winds  or  radiation;  or 
are  accreted  onto  larger  objects  or  the  central  star  itself.  Observational  evidence  suggests 
that  most  disks  effectively  disappear  by  about  10  Myr  (Beckwith,  1999). 

Many  physical  processes  are  involved  in  the  transformation  of  composition  and  struc¬ 
ture.  This  chapter  first  discusses  how  circumstellar  disks  are  formed  during  the  star  forma¬ 
tion  process  and  then  discusses  the  processes  thought  to  play  the  greatest  roles  in  deter¬ 
mining  the  structure  and  content  of  disks. 

2.2  Formation  of  Circumstellar  Disks 

The  current  paradigm  of  star  formation  (described  by  Shu,  Adams,  and  Lizano  1987) 
describes  a  scenario  in  which  stars  form  from  the  collapse  of  slowly  rotating  molecular  cloud 
cores.  The  onset  of  the  collapse  is  a  result  of  a  density  perturbation  within  a  cloud  core. 
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Perturbations  are  generally  thought  to  be  the  result  of  inhomogeneities  within  the  cloud  or 
the  result  of  an  outside  disturbance.  Recent  evidence  from  analysis  of  meteorites  suggests 
that  the  solar  system  was  created  when  a  cloud  was  perturbed  by  a  nearby  supernova 
(Yulsman  2003).  The  collapse  of  a  core  leads  to  a  protostar  and  a  flattened  disk  surrounded 
by  an  envelope  of  dust  and  gas.  Theory  predicts  that  the  disk  forms  as  a  result  of  the  initial 
rotation  of  the  cloud  core  and  the  conservation  of  angular  momentum  during  the  collapse. 
The  material  in  the  envelope  is  generally  removed  from  the  system  through  an  amalgam  of 
processes.  These  clearing  processes  include  accretion  onto  the  star,  stellar  winds,  radiation 
pressure,  and  accretion  onto  growing  planetesimals.  The  planetesimals  that  are  initially 
created  by  mutual  collisions  of  dust  grains  eventually  lead  to  the  large  planet  sized  bodies 
that  populate  a  planetary  system. 

The  following  sections  discuss  in  more  detail  the  processes  that  lead  to  the  formation 
of  the  disk  and  its  structure;  the  formation  of  kilometer  sized  bodies  from  micron  sized  dust 
grains;  and  Anally  the  formation  of  terrestrial  and  giant  planets. 

2.2.1  Formation  and  Structure  of  a  Protoplanetary  Disk 

The  star  formation  process  begins  in  the  dense  regions  of  giant  molecular  clouds 
(GMCs).  In  our  galaxy,  GMCs  exist  throughout  the  galaxy  but  are  primarily  found  in  the 
spiral  arms.  They  are  cold  clouds  of  primarily  hydrogen  gas  that  have  masses  10^  M0 
and  are  gravitationally  bound  (Blitz  and  Williams  1999).  Within  the  GMCs  are  smaller, 
but  still  gravitationally  bound,  regions  known  as  cores.  It  is  from  these  cores  that  the 
individual  stars  are  likely  to  form.  These  cores  have  typical  masses  of  ~  10  Mq  and  typical 
sizes  of  0.1  pc  (Shu  et  al.  1999).  They  are  generally  supported  through  a  combination  of 
magnetic  flelds,  turbulence,  and  gas  pressure.  When  the  mass  of  a  core  reaches  a  critical 
value,  known  as  the  Jean’s  mass,  the  core  begins  to  contract  and  the  star  formation  process 
begins.  The  Jean’s  mass  is  estimated  by  considering  kinetic  and  potential  energy.  The 
magnetic  flelds,  turbulence,  and  gas  pressure  also  play  a  role  in  determining  the  critical 
value. 

As  the  irregularly  shaped  core  of  dust  and  gas  collapses  and  begins  forming  a  dense 
central  object,  the  system  must  account  for  the  angular  momentum  associated  with  the 
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original  rotating  core.  Some  material  is  accreted  onto  the  central  protostar  but  other 
material  retains  enough  of  its  original  angular  momentum  to  settle  into  stable  orbits  and 
create  a  thin  protostellar  disk.  The  disk  is  expected  to  form  roughly  10^  yrs  after  the 
beginning  of  the  initial  freefall  collapse  (Beckwith  1999). 

A  protoplanetary  disk  is  an  example  of  an  accretion  disk  which  is  defined  as  a  flattened 
disk  composed  of  material  that  is  spiraling  down  onto  the  surface  of  a  central  mass.  Through 
the  conservation  of  angular  momentum,  the  material  within  an  accretion  disk  is  transported 
inward  while  the  angular  momentum  is  transferred  outward.  This  transport  occurs  as  a 
result  of  the  viscous  efiects  of  the  material  within  the  disk  (Kenyon  1999).  As  frictional 
forces  cause  heating  of  material  due  to  differential  rotation,  some  material  moves  due  to 
net  loss  of  kinetic  energy.  This  motion  also  results  in  a  decrease  in  angular  momentum 
and  this  must  be  balanced  by  outward  motion  of  other  material  to  conserve  total  angular 
momentum  within  the  disk. 

The  transport  of  material  inward  and  angular  momentum  outward  plays  a  key  role  in 
defining  the  structure  and  evolution  of  the  protoplanetary  disk.  As  mass  is  moved  inward 
and  accreted  on  to  the  star  or  possibly  removed  by  outflows,  the  surface  density  spreads 
outward  (Ruden  1999) .  This  results  in  a  dispersal  of  the  disk  material  over  time. 

Since  it  is  not  possible  to  directly  image  most  disks,  the  general  structure  of  disks 
has  been  gleaned  through  analysis  of  the  spectral  energy  distributions  (SEDs)  of  star-disk 
systems.  SEDs  show  the  measured  flux  from  a  somce  as  a  function  of  wavelength.  SEDs 
associated  with  a  star  surrounded  by  a  disk  show  particular  characteristics  that  clearly 
confirm  both  the  presence  and  the  details  of  the  disk.  The  SED  for  a  star  is  generally  that 
of  a  black  body  at  the  effective  temperature  of  the  star.  The  SED  for  a  star  and  a  disk 
will  show  not  only  the  black  body  curve  for  the  star,  but  also  an  excess  flux  in  the  infrared 
region  of  the  SED.  This  excess  is  associated  with  thermal  emissions  firom  the  dust  in  the 
disk  that  is  heated  primarily  by  radiation  from  the  central  star. 

As  the  star  proceeds  through  its  early  evolutionary  stages,  the  SEDs  also  evolve.  In 
the  early  stages  the  SED  is  actually  dominated  by  the  original  envelope  of  dust  and  gas. 
As  the  star  evolves  material  is  accreted,  is  blown  out  by  radiation  pressure,  or  settles  into  a 
circumstellar  disk.  As  this  occurs,  the  SED  changes  to  one  in  which  the  black  body  curve 
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of  the  star  can  be  seen  but  is  modified  by  a  superposition  of  black  body  curves  associated 
with  a  range  of  dust  particles  at  lower  temperatures. 

Mass  estimates  of  disks  are  most  commonly  made  through  long  wavelength  measure¬ 
ments  (A  >  300tJ.m).  Typical  masses  for  disks  determined  from  sub-mm  observations  are 
on  the  order  of  a  few  percent  of  a  solar  mass  (Beckwith  1999).  This  fact  is  significant  as 
it  is  considerably  higher  than  the  mass  thought  necessary  to  form  all  of  the  planets  within 
our  solar  system. 


2.2.2  Formation  of  Planetesimals 

The  initial  protoplanetary  disk  is  almost  entirely  (~99%)  hydrogen  and  helium  gas. 
However,  as  the  protostar  evolves,  the  gas  is  removed  from  the  disk  primarily  by  radiation 
pressure  and  stellar  winds.  During  this  time  the  dust  grains,  that  make  up  the  remaining 
~1%  of  the  disk,  begin  the  process  of  growth  to  larger  planetesimals. 

While  gas  is  still  present  within  the  disk,  solid  particles  are  affected  by  gas  drag.  The 
effects  of  gas  drag  are  dependent  on  the  size  of  the  particles.  The  smaller  particles,  which  are 
coupled  to  the  gas,  move  at  sub-Keplerian  velocities  while  the  larger  particles,  which  are  less 
affected  by  the  gas  drag  force,  continue  to  move  at  near  Keplerian  velocities  but  experience 
orbital  decay  (Ruden  1999).  This  results  in  a  differential  inward  drift  that  increases  the 
firequency  of  collisions  between  particles.  Many  of  these  particles  stick  together  forming 
larger  particles.  This  process  marks  the  beginning  of  the  growth  to  larger  objects  within 
the  disk. 

As  particles  grow  through  collisions  induced  by  the  inward  drift  they  also  begin  settling 
into  the  mid-plane  of  the  disk.  Settling  occurs  because  the  mid-plane  defines  the  minimum 
in  the  vertical  gravitational  potential.  This  is  described  through  the  vertical  acceleration 
measured  with  respect  to  the  disk,  and  is  expressed  as 
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(2.1) 


where  z  measures  the  vertical  distance  from  the  mid-plane  and  Q  is  the  orbital  frequency 
at  a  given  radius  from  the  star.  The  above  expression  implies  that  the  smaller  grains,  with 
a  longer  orbital  period  due  to  the  gas  drag,  will  take  the  longest  to  settle  (Ruden  1999). 
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This  leads  to  a  picture  in  which  the  larger  particles  will  settle  to  the  disk  mid-plane  first 
and  begin  the  process  of  further  growth  as  the  mid-plane  density  increases.  The  process 
is  enhanced  as  the  particles  get  larger  and  present  a  larger  collisional  cross-sectional  area. 
The  final  result  of  this  process  is  a  thin  layer  of  large  particles  in  the  disk’s  mid-plane. 

As  the  particles  become  larger,  they  become  less  affected  by  gas  drag  and  eventually 
reach  Keplerian  orbits.  This  occurs  when  the  objects  have  grown  from  small  dust  particles 
to  bodies  with  diameters  ^  1  km,  generally  referred  to  as  planetesimals.  As  the  small 
particles  continue  to  drift  inward  past  the  larger  planetesimals,  they  are  accreted  on  to 
these  large  bodies  and  are  thought  to  be  capable  of  producing  objects  on  the  scale  of  ~  100 
km  on  timescales  of  ^  10®  years  (Ruden  1999). 

It  is  from  this  distribution  of  planetesimals  scattered  throughout  the  disk  that  planets 
begin  to  grow  through  the  accretion  process. 

2.2.3  Formation  of  Planets 

Planetesimals  distributed  through  the  mid-plane  of  a  disk  continue  to  experience  colli¬ 
sions  and  create  larger  bodies.  As  the  bodies  become  more  massive  the  effective  collisional 
cross  section  increases  as  gravitational  effects  begin  to  become  important.  This  effect  can 
be  described  by  the  gravitational  focusing  factor  expressed  as  1  -b  where  V  is  the 

relative  collisional  velocity.  The  term  Ug  is  the  mutual  escape  velocity  and  is  expressed  as 

Ve  = 

where  mi  and  m2  are  the  masses  of  the  two  bodies  and  ai  and  02  are  the  radii.  Through  this 
gravitational  focusing,  bodies  continue  to  grow  and  simultaneously  increase  their  collisional 
cross  section  thus  enabling  them  to  accelerate  their  growth. 

The  idea  of  runaway  growth  has  been  considered  as  a  mechanism  which  can  increase 
the  rate  of  growth  for  a  protoplanet  during  the  process  of  planet  formation.  Runaway 
growth  is  the  result  of  an  increase  in  the  gravitational  focusing  factor  which  occurs  because 
of  the  increase  in  mass.  This  occurs  when  the  ratio  of  /V"^  becomes  very  large,  ^  1. 


12 


This  critical  ratio  is  driven  both  by  the  size  of  the  planetesimal  and  by  the  relative 
velocity  V.  The  relative  velocity  is  influenced  primarily  by  gravitational  scattering,  damping 
due  to  gas  drag,  and  energy  transfer  via  dynamical  friction  (Ruden  1999) .  A  planetesimal 
will  continue  to  grow  through  this  process  until  it  has  accumulated  all  material  in  what  is 
called  the  accretion  zone.  The  accretion  zone  is  deflned  by  the  Hill  sphere  which  is  based 
on  determining  the  equilibrium  point  where  tidal  effects  from  the  star  are  balanced  by  the 
gravity  of  the  planet.  The  radius  of  the  Hill  sphere  is  given  as 


where  r  is  the  distance  from  the  star,  rrip  is  the  mass  of  the  planet,  and  M*  is  the  mass 
of  the  star.  It  is  generally  thought  that  the  accretion  zone  is  approximately  four  times  the 
radius  of  the  Hill  sphere  (Lissauer  1993). 

If  a  reservoir  of  gas  exists  near  the  forming  planet,  the  gas  will  be  accreted.  Accretion 
will  occur  provided  a  critical  mass  of  the  core  is  achieved.  This  critical  mass  is  determined 
by  finding  the  mass  at  which  gravity  precludes  hydrostatic  equilibrium.  This  critical  mass 
is  generally  thought  to  be  ;^10-15M©  (Ruden  1999).  The  accretion  of  gas  onto  these  rocky, 
critical  mass  cores  is  the  process  that  is  thought  to  lead  to  the  formation  of  the  gas  giants. 

Other  process  also  occur  during  the  planet  formation  process  that  contribute  to  the 
structure  of  the  planetary  system.  These  include  tidal  interactions  between  the  disk  and 
the  planets  that  lead  to  both  gap  formation  and  planetary  migration.  This  will  be  discussed 
later  in  chapters  4  and  8. 

2.3  Physical  Processes  Within  a  Circumstellar  Disk 

The  primary  physical  processes  which  occur  within  a  circumstellar  disk  that  influence 
the  dynamics  of  the  constituent  particles  are  associated  with  gravitational  perturbations, 
radiation,  interaction  with  gas,  and  collisions.  With  each  of  these  processes  is  an  associated 
time  scale  that  is  based  on  environmental  conditions  influencing  the  particular  process.  An 
analysis  of  the  timescales  is  important  for  two  reasons.  First,  for  processes  that  occur  on 
long  time  scales,  it  is  necessary  to  verify  that  the  system  has  had  enough  time  to  fully 
evolve  under  the  influences  of  that  particular  process.  Second,  the  timescales  need  to  be 
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compared  to  one  another  in  order  to  see  which  process  will  dominate  over  others  within  the 
disk.  In  general,  the  processes  that  occur  on  the  shortest  timescales  will  have  the  greatest 
effect  on  the  disk.  For  example,  if  a  certain  sized  particle  is  quickly  destroyed  by  collisions, 
such  as  the  case  for  micron  sized  particles  in  Beta  Pictoris,  the  dust  particles  do  not  have 
time  to  spiral  into  the  star.  This  implies  that  the  collisional  processes  will  dominate  and 
play  a  significant  role  in  influencing  the  distribution  of  these  particular  particles. 

This  section  will  give  a  brief  description  of  each  of  the  important  processes  that  occur 
within  a  circumstellar  disk  and  discuss  their  corresponding  timescales. 


2.3.1  Gravity 

The  dominant  influence  in  the  circumstellar  environment  is  that  of  gravity  associated 
with  the  central  star.  Gravitational  perturbations  may  also  come  from  other  planets  within 
the  system.  These  perturbations  resulting  from  the  secondary  masses  are  what  drive  the 
secular  effects  and  will  be  discussed  in  detail  in  the  next  section.  The  timescale  associated 
with  the  influences  of  the  central  star  is  simply  the  Keplerian  period.  This  is  determined 
from  Kepler’s  third  law  that  states  the  square  of  the  orbital  period  is  proportional  to  the 
semi-major  axis  cubed.  This  law  is  generally  expressed  as 


(2.4) 


where  p  =  G{mi  +m2).  This  can  also  be  expressed  as  a  mean  rate  in  terms  of  solar  system 
parameters,  i.e.  AUs  and  years. 
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where  tUc  and  m  are  in  units  of  solar  masses,  a  is  in  AU,  and  n  is  in  units  of  radians  per 
year.  This  can  also  be  expressed  as  a  characteristic  timescale 


iper  — 


(2.6) 


where  a  is  again  in  AUs  and  rric  and  m  are  in  units  of  solar  masses. 
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These  values  of  frequency  and  time  provide  a  standard  of  measure  used  in  evaluating 
the  frequency  and  timescales  of  the  other  processes  described  below. 


2.3.2  Radiative  Forces 

All  particles  within  a  circumstellar  disk  interact  with  radiation  from  the  central  star. 
The  magnitude  of  the  radiation  effects  is  based  on  the  characteristic  cross  sectional  area  of 
a  particle  and  its  mass.  For  larger  particles,  gravity  dominates  over  the  radiative  forces. 
However,  for  smaller  particles  the  ratio  of  radiative  forces  to  gravitational  forces  approaches 
or  exceeds  unity.  This  ratio  can  be  expressed  with  the  dimensionless  quantity  p  (see  for 
example  Gustafson  1994  and  references  within).  The  value  of  ^  can  be  expressed  as 
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(2.7) 


where  So  is  the  stellar  flux  density  and  c  is  the  speed  of  light.  Parameterizing  ^  in  terms  of 
solar  values  uses  the  constant  Cr  =  7.6  x  10^  g  cm?  and  defines  cr/m  as  the  cross-sectional- 
area-to-mass  ratio  for  a  particle  of  diameter  D  and  density  p.  The  quantity  Qpr  is  the 
radiation  pressure  efficiency  and  is  related  to  the  absorption  and  scattering  eflBiciencies; 
Qpr  =  Qabs  +  Qaca-  For  dust  particles  assumed  to  be  homogeneous  spheres,  Qpr  is  usually 
calculated  from  Mie  theory  (Gustafson  1994).  Sample  calculations  of  ^8  as  a  function  of 
particle  diameter  are  shown  in  figme  2-1.  The  optical  properties  and  densities  axe  those 
of  the  particles  considered  to  most  likely  be  representative  of  the  dust  surrounding  Beta 
Pictoris.  This  will  be  discussed  further  in  chapter  7. 

Radiation  is  directed  radially  away  from  the  star  and  therefore  radiation  acts  to  cormter 
the  effect  of  gravity.  This  force  effectively  reduces  the  mass  of  the  star  from  the  perspective 
of  the  dust  particle.  Equation  2.7  shows  that  is  dependent  on  the  size  of  the  particle 
through  the  a/m  term.  This  imphes  an  inverse  relationship  between  S  and  the  particle 
diameter,  i.e.  P  oc  1/D.  Since  /3  is  size  dependent,  particles  of  different  diameter  Z),  will  see 
a  different  potential.  This  means  that  particles  with  the  same  initial  conditions  of  position 
and  velocity,  but  with  different  diameters,  will  move  on  different  orbits. 
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Particle  Size  (microns) 


Figure  2-1:  The  parameter  Beta  in  Beta  Pictoris  (M  =  l.SM©;  L  =  S.TL©)  as  a  function 
of  particle  size,  calculated  using  values  for  Qpr  determined  from  Mie  theory.  Two  different 
classes  of  materials  are  shown:  materials  described  as  astronomical  silicates  (Laor  and 
Draine,  1993)  and  organic  refractory  material  Li  and  Greenberg  (1997).  For  this  research, 
these  materials  are  the  primary  candidates  for  simulating  the  dust  within  Beta  Pictoris. 
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Particles  can  thus  be  grouped  into  different  subcategories  based  on  their  P  value.  Lcirge 
particles  have  a  /0  value  that  approaches  zero  and  therefore  are  relatively  unaffected  by 
radiation  pressure.  The  smallest  particles,  the  0  meteoroids,  have  the  largest  0  values, 
greater  than  0.5,  and  are  quickly  expelled  from  the  system  on  hyperbolic  orbits.  Particles 
that  have  intermediate  values,  0.0  <  <  0.5,  exist  on  bound  but  eccentric  orbits  and  are 

referred  to  as  a  meteoroids. 

Particles  feel  another  effect  of  radiation  resulting  from  the  component  of  the  radiation 
force  tangential  to  the  particles  orbit.  This  force  acts  as  a  dynamic  drag  force  and  is  referred 
to  as  Poynting-Robertson  (PR)  drag.  PR  drag  acts  to  decrease  both  the  eccentricity  and 
the  semi-major  axis  in  a  manner  proportional  to  0  (Burns  et  al.  1979).  This  drag  force 
causes  the  particle,  initially  at  ri,  to  spiral  toward  the  star  and  reach  a  distance  r2  on 
timescales  that  are  a  function  of  0  (Wyatt  et  al.  1999): 


where  tpR  is  in  years.  It  is  important  to  note  that  tpR  is  inversely  proportional  to  0  and 
therefore  proportional  to  the  particle  size,  D,  as  this  effect  further  differentiates  the  behav¬ 
ior  of  particles  based  on  size. 

2.3.3  Collisions 

Collisions  play  important  roles  in  the  circumstellar  environment  as  they  provide  mech¬ 
anisms  for  both  creation  and  destruction  of  dust  particles.  Based  on  collisional  models  of 
dust  in  disk  environments,  the  size  distribution  of  particles  is  rather  steep,  i.e.  there  are 
many  more  small  particles  than  large  particles.  Due  to  the  steep  size  distribution  expected 
for  particles  in  the  disk  environment,  the  surface  area  of  the  disk  is  dominated  by  the  smaller 
particles  (for  a  thorough  discussion  see  Wyatt  et  al.  (1999)).  These  particles,  however,  may 
have  lifetimes  associated  with  radiation  effects  that  are  much  shorter  than  the  age  of  the 
disk,  i.e.  these  particles  cannot  be  left  over  from  the  original  formation  of  the  disk.  This 
condition  necessitates  a  source  through  which  these  particles  can  be  replenished. 

The  timescale  associated  with  the  destruction  of  the  dust  particles  is  based  on  the 
assumption  that  a  typical  particle  size  of  the  emitting  dust  can  be  determined.  The  typical 


17 


size  is  based  on  assumptions  about  the  particles  responsible  for  the  majority  of  the  cross 
sectional  area  observed  at  a  particular  wavelength.  This  is  estimated  by  considering  both 
number  density  and  the  physical  cross-section  of  the  particle.  In  terms  of  this  typical 
particle  size,  Wyatt  et  al.  (1999)  provided  an  expression  for  the  collisional  timescale  as 

^colliDtypjf)  —  tper  (r)/47rre//(r)  (2.9) 

where  Dtyp  is  the  characteristic  size  of  the  particles  responsible  for  most  of  the  disk’s  cross 
sectional  area.  The  optical  depth  can  be  thought  of  as  a  measure  of  the  total  material  in 
a  particular  line  of  sight  modulo  the  effects  of  absorption  and  scattering  properties  of  the 
dust.  These  optical  properties  are  quantified  by  the  extinction  coefficient.  The  effective 
optical  depth  in  the  above  expression  is  the  face-on  optical  depth  of  a  disk  observed  if  the 
extinction  coefficients  for  the  particles  are  equal  to  1.  The  face-on  optical  depth  can  be 
estimated  observationally  by  considering  the  simplified  relationship  (Wyatt  et  al.  1999) 

J7j/(A,r)  ~  Tg^^(r)floi,s5i/(A,  iAtyp,  r)Qa6s('^) -^tj/p)  (2.10) 

which  relates  the  flux  density  received  at  earth,  F^,  at  a.  given  wavelength.  A,  from  a  particle 
with  a  Planck  function  of  at  a  distance  r  from  the  star,  and  with  characteristic  size 
Dtyp.  The  observational  solid  angle  is  related  to  the  detector  pixel  size,  Q  ~  dpiJ^,  and  the 
absorption  coefficient  of  the  typical  particle  is  Qabs-  Since  disks  are  unlikely  to  be  observed 
directly  face-on,  projection  effects  must  be  considered  to  transform  Tg/f  into  a  face-on  Tgff. 
Once  the  face-on  or  normal  optical  depth  is  determined,  the  collisional  timescales  can  be 
estimated. 

2.3.4  Other  Forces 

The  circumstellar  environment  is  very  complex  and  has  many  other  factors,  besides  the 
ones  already  mentioned,  that  influence  the  dynamics  of  the  dust  particles  within  the  disk. 
The  effects  that  often  play  significant  roles  in  this  environment  are  briefly  described  below 
along  with  the  justification  for  not  emphasizing  their  impact  on  the  modeling  performed  in 
this  research. 
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Stellar  wind  effects 

The  stellar  wind  will  interact  with  the  dust  particles  in  two  primary  ways.  First, 
because  the  dust  particles  are  charged  and  because  the  winds  carry  with  them  a  magnetic 
field,  a  Lorentz  force  arises.  The  Lorentz  force  is  a  function  of  the  overall  charge  of  the 
dust  particle  and  the  strength  of  the  interacting  magnetic  field.  Second,  actual  collisions 
between  the  dust  particles  and  the  protons  comprising  the  wind  give  rise  to  an  additional 
drag  component  (see  Gustafson  (1994)  for  example). 

Stellar  winds  are  thought  to  be  associated  with  coronal  holes  and  are  thus  tied  to 
coronal  and  possibly  chromospheric  activity.  A-tjrpe  stars  do  not  possess  a  deep  convective 
zone  from  which  high  speed  winds  are  effectively  driven  (Panzera  et  al.,  1999).  Since  Beta 
Pictoris  is  an  A-type  star  and  since  an  A-type  star  is  expected  to  have  little  chromospheric 
and  coronal  activity  (Artymowicz  1997)  the  effects  of  the  stellar  winds  around  Beta  Pictoris 
are  neglected  in  this  research. 

Gas  drag 

For  very  young  disks  in  which  the  gas  to  dust  mass  ratio  is  significant,  the  dynamic 
drag  on  the  dust  particles  caused  by  the  gas  must  be  considered.  However,  Liseau  and  Arty¬ 
mowicz  (1997)  determined  from  mm-wavelength  observations  that  the  overall  gas  to  dust 
mass  ratio  was  <  0.003  and  that  this  indicates  the  disk  is  effectively  gas  free  (Artymowicz, 
1999).  Therefore,  gas  drag  within  the  Beta  Pictoris  system  was  assumed  to  be  negligible 
and  this  effect  was  not  considered  any  further  in  this  research. 

2.3.5  Conclusions  on  Timescales 

The  overall  picture  that  emerges  from  this  analysis  for  an  evolved  circumstellar  disk  is 
that  the  smallest  particles  are  quickly  removed  from  the  system  through  radiation  pressure. 
The  medium  sized  particles  are  either  removed  by  collisions  or  spiral  into  the  star  through 
the  effects  of  PR  drag.  The  largest  particles,  which  axe  responsible  for  the  creation  of 
the  smaller  particles  through  collisional  cascade,  are  affected  by  gravitational  forces  alone. 
Therefore  it  can  be  concluded  that,  in  general,  the  distribution  of  larger  bodies  within  a 
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circumstellar  disk  is  primarily  determined  by  the  sum  of  all  gravitational  effects  within  the 
disk.  These  gravitational  effects  can  result  from  the  central  star,  other  large  bodies,  and 
possibly  from  the  disk  itself  provided  it  has  a  significant  mass  compared  to  the  other  masses 
within  the  system.  This  picture  will  be  discussed  further  in  chapter  7  as  it  applies  to  Beta 
Pictoris. 

The  following  chapter  discusses  the  method  of  determining  the  combined  effects  of 
gravity  in  a  complex  planetary  system  such  as  a  circumstellar  disk  with  multiple  imbedded 
planets.  This  will  be  done  through  the  theoretical  framework  of  secular  perturbation  theory. 


CHAPTER  3 

SECULAR  PERTURBATION  THEORY 


3.1  Introduction 

The  previous  chapter  began  to  reveal  the  many  influences  that  play  a  part  in  deter¬ 
mining  the  characteristics  of  a  circumstellar  disk.  These  influences  were  associated  with 
such  things  as  radiation-particle  interaction,  collisions,  drag  forces,  etc.  This  chapter  will 
expand  on  the  discussion  by  providing  a  summary  of  the  theory  used  to  describe  the  effects 
of  gravity  from  not  only  the  central  star,  but  from  other  planets  in  the  system. 

Generally  speaking,  a  circumstellar  disk  is  a  flattened  cloud  of  dust  in  which  each  dust 
particle  moves  on  an  its  own  individual  orbit.  In  order  to  understand  the  effects  of  planets 
on  the  disk,  a  link  must  be  established  between  the  planets  and  the  dust  particle  orbits. 
The  theoretical  framework  through  which  this  is  accomplished  in  this  research  is  secular 
perturbation  theory.  Secular  perturbation  theory  is  a  well  established  approach  that  has 
been  applied  extensively  to  problems  within  the  solar  system.  In  fact,  analysis  of  the  secular 
perturbations  arising  from  Jupiter  and  Satmrn  can  be  shown  to  create  a  forced  inclination 
that  varies  with  semi-major  axis,  a  -  see,  for  example,  Dermott  et  al.  (1984,  1985)  and 
Murray  and  Dermott  (1999).  This  is  significant  as  this  characteristic  is  the  source  of  the 
mid-plane  warping  that  is  the  crux  of  this  research. 

This  chapter  will  explain  how  expressions  describing  particle  orbits  can  be  developed  by 
considering  the  gravitational  perturbations  from  massive  bodies.  The  solutions  that  result 
from  this  approach  describe  the  time  evolution  of  each  of  the  particles  orbital  elements: 
eccentricity,  longitude  of  pericenter,  inclination,  and  longitude  of  ascending  node.  A  very 
important  concept  that  will  be  repeated  throughout  this  discussion  is  that  the  solutions 
can  be  broken  down  into  two  components.  These  components  are  referred  to  as  the  proper 
and  forced  orbital  elements.  The  proper  elements  are  related  to  the  initial  conditions  at 
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creation  of  the  particle  under  consideration.  The  forced  elements  correspond  to  the  effects 
of  the  perturbations  and  are  thus  linked  to  the  planets. 

The  goal  of  this  chapter  is  to  explain  the  theoretical  framework  that  produces  the 
forced  orbital  elements  from  a  particular  planetary  configuration.  It  is  the  forced  orbital 
elements  that  are  then  used  as  an  input  to  the  modeling  process  described  in  the  later 
chapters. 


3.2  Second  Order  Secular  Perturbation  Theory 
Since  the  many-body  problem  in  dynamics  is  insoluble  analytically,  alternatives  to  an 
exact  solution  have  been  developed.  The  approach  outlined  here  is  to  start  with  the  three- 
body  problem  and  then  divide  the  resulting  potential  into  two  components,  one  of  which 
represents  the  two-body  central  mass  and  the  second  of  which  represents  the  portion  of  the 
potential  arising  from  the  additional  masses.  This  second  component  is  known  as  the  Dis¬ 
turbing  Function.  A  detailed  derivation  of  the  disturbing  function  is  presented  in  Murray 
and  Dermott  (1999). 

3.2.1  The  Disturbing  Function 

Starting  with  the  three-body  equations  of  motion  in  cartesian  coordinates,  the  accel¬ 
eration  of  the  secondary  masses  with  respect  to  the  central  mass  can  be  expressed  as 

Ti  =  ViiUi+Ui)  (3.1) 

rj  =  ^ j{Uj -i-TZj) 

where  the  potential  has  been  broken  up  into  two  components;  one  describing  the  two  body 
part,  U,  where 

a  (”»+”“)  (3.2) 

n 

Q  (tHc  -k  mj) 


Ui  = 

= 


and  one  describing  the  portion  of  the  potential  arising  from  the  gravitational  perturbations 
fi-om  the  secondary  masses,  TZ.  This  portion  of  the  potential  is  referred  to  as  the  Disturbing 
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Function.  For  the  three  body  problem,  the  disturbing  function  for  the  individual  masses 
can  be  expressed  as 


Ui  =  G 


Kj  =  G 


Tj  -  Ti 
mj 

Ti-Tj 


(3.3) 


Gmi 

rf 


In  order  to  handle  the  rj  —  rj  terms  analytically,  the  disturbing  function  is  usually 
expanded  in  terms  of  Legendre  polynomials.  Furthermore,  the  disturbing  function  can  be 
transformed  from  a  form  based  on  positions  and  velocities  to  a  form  expressed  in  terms 
of  the  osculating  orbital  elements,  a,  e,  I,  A,  zo,  and  fi  which  represent  the  semi-major 
axis,  eccentricity,  inclination,  mean  longitude,  longitude  of  pericenter,  and  longitude  of  the 
ascending  node  respectively,  i.e. 


TZ  =  f{m,  m',  a,  a',  e,  e',  I,  I’,  cos  <f>) 


where  the  argument  of  the  cosine  function  is  a  linear  combination  of  the  mean  longitudes, 
longitude  of  pericenter,  and  longitude  of  ascending  node.  This  term  is  expressed  as 

<f>  =  A'  +  jiX  +  hzu'  +  jiw  +  j^O!  +  j&Q.  (3.4) 


where  the  unprimed  quantities  represent  the  inner  mass  and  the  primed  quantities  represent 
the  outer  mass. 


3.2.2  The  Secular  Portion  of  the  Disturbing  Function 

The  terms  within  the  disturbing  function  can  be  grouped  together  by  characteristic 
timescales.  The  three  categories  are  short  period  terms,  mean  motion  resonance  terms, 
and  secular  terms.  The  short  period  terms  axe  dependent  on  the  mean  longitude  A,  and 
represent  effects  that  occur  on  the  scale  of  orbital  periods.  The  terms  associated  with  mean 
motion  resonances  represent  effects  that  axe  also  dependent  on  the  mean  longitude  and 
would  only  effect  certain  locations  within  a  disk,  generally  at  radii  where  the  corresponding 
period  has  a  commensurability  with  one  of  the  planets.  The  secular  terms  produce  effects 
that  are  felt  throughout  the  disk  and  are  independent  of  mean  longitude.  The  independence 
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of  mean  longitude  implies  that  the  instantaneous  positions  of  the  masses  have  no  effect  on 
large  timescales.  For  this  research,  long  term  effects  that  are  felt  throughout  the  disk  will 
be  considered  as  the  source  of  perturbations  leading  to  large  scale  warping,  and  thus  only 
the  secular  part  of  the  disturbing  function  will  be  retained. 

The  time  variation  of  the  orbital  elements  can  be  expressed  through  the  use  of  La¬ 
grange’s  Planetary  Equations: 


ej 

1  dUj 
nja^ej  dwj 

(3.5a) 

Wj 

1  mj 

UjO.'jej  dej 

(3.5b) 

ij 

1  dUj 
njapj  dQj 

(3.5c) 

tlj 

1 

II 

(3.5d) 

This  particular  form  of  Lagrange’s  equations  assumes  low  inclination  and  eccentricity. 

This  research  will  consider  massless  test  particles  disturbed  by  an 

arbitrary  number 

of  planets  and  the  resulting  effects  on 

the  large  scale  structure  of  a  disk 

.  In  this  case,  the 

disturbing  function  can  be  expressed  as 

TZ  =  + 

-BI^ 

2 

(3.6) 

N 

N 

+  AjCCj  cos 

(tx7  —  Wj)  +  Bjllj  COs(J7  —  fij)  1 

j=l 

where  the  values  of  A  and  B  are  functions  of  mass  and  semi-major  axis  only,  and  are  given 

by 

A  =  - 

(3.7a) 

Aj  = 

(3.7b) 

B  =  - 

i=i 

(3.7c) 

Bj  =  - 

1  ^  ,  (1)  /  \ 

(3.7d) 
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The  value  of  a  is  determined  through  the  ratios  of  semi-major  axes  of  the  test  particle  and 
the  perturbing  mass.  More  precisely  aj  =  aj/a  if  aj  <  a  and  a/aj  if  aj  >  a.  Similarly, 
ocj  =  aj/a  if  aj  <  a  and  a/cj  if  aj  >  a.  The  quantities  b3y2(«)  axe  the  Laplace 

coefficients  and  are  given  by  the  following  expressions 


1  f  cos  ipdtp 

J  (1  —  2q:cosV’ +  Q!^)5 

27r 

1  r  cos  2'tl)d‘tp 

(1  —  2q:cosV’ + 


(3.8a) 


(3.8b) 


The  form  of  the  disturbing  function  presented  here  is  to  second  order  in  inclination 
and  eccentricity  eej,  and  II j.  The  truncation  of  the  disturbing  function  limits  the 

analysis  to  small  eccentricities  and  inclinations. 


3.3  Solutions  to  Lagrange’s  Equations 

The  form  of  Lagrange’s  planetary  equations  shown  in  the  previous  section  has  obvious 
difficulties  when  the  eccentricity,  e,  and  inclination,  I,  are  small.  To  remove  the  singular¬ 
ities  when  e  and  I  are  small,  the  different  components  can  be  used  together  in  a  vector 


representation  with  the  following  components 

h  =  esmw  (3.9a) 

k  =  ecosm  (3.9b) 

p  =  /sinfi  (3.9c) 

q  =  /cosfi.  (3.9d) 


These  variables  define  the  components  of  the  eccentricity  and  inclination  vectors  and  are 
functions  of  the  osculating  orbital  parameters  -  the  instantaneous  parameters  determined 
for  a  specific  particle  trajectory.  Using  these  components  of  the  eccentricity  and  inclination 
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vectors  the  disturbing  function  takes  the  form 

n  =  na^^A{h’^  +  k^)  +  ^B{p'^  +  q^)  (3.10) 

N  N 

+  Aj {hhj  +  kkj)  +  Y  ^3  (PPj  +  } 

j=i  j=i 

Likewise,  Lagrange’s  planetary  equations  take  on  the  form 


h  = 

1  an 

TijOj  dk 

(3.11a) 

k  = 

1  an 

rijOj  ah 

(3.11b) 

P  = 

1  an 

rijoj  aq 

(3.11c) 

Q  = 

1  an 

njaj  ap 

(3.11d) 

Solving  these  equations  provides  solutions  to  the  equations  of  motion  in  terms  of  the 
components  h,  k,  p,  and  q.  The  solutions  to  these  equations  are 


h  =  Cproper  sin(^t  +  ^)  +  ho{t)  (3.12a) 

k  =  eproper  cos{At  +  ^)  +  ko{t)  (3.12b) 

P  —  Iproper  sin(5t  +  7)  +po(<)  (3.12c) 

q  =  Iproper  cos{Bt  +  7)  +  qo{t)  (3.12d) 


where  Cproper  and  Iproper  are  the  proper  eccentricity  and  proper  inchnation,  or  the  eccentric¬ 
ity  and  inclination  associated  with  the  initial  conditions  of  the  test  particles.  The  quantities 
At  +  P  and  Bt  +  'y  are  equal  to  Wproper  and  ^proper  •>  the  test  particles  proper  pericenter  and 
longitude  of  ascending  node  respectively.  As  was  stated  earlier  in  the  discussion,  the  two 
components  of  the  solutions  eure  related  to  the  two  primary  aspects  of  the  problem.  The 
first  term  is  related  to  the  initial  conditions  of  the  massless  particle  under  consideration. 
The  second  term  in  each  equation  is  associated  with  the  perturbation  and  can  be  expressed 


26 


as 


o 

II 

X]  A  sin{git  +  Pi) 

(3.13a) 

ko{t)  = 

N 

XI  A  cos{git  +  Pi) 

(3.13b) 

Po{t)  = 

H  f-  +  T'*) 

i=\ 

(3.13c) 

II 

Xd^  r.  COs(/it  +  7i) 

i=l  ^  •'* 

(3.13d) 

where  gi  and  fi  are  system  eigenfrequencies.  Eigenfrequencies  are  characteristic  frequencies 
associated  with  all  dynamical  systems  possessing  internal  periodicity.  It  is  the  eigenfrequen¬ 
cies  of  a  system  that  determine  at  which  frequencies  the  secular  resonances  will  occur.  It 
is  important  to  remember  that  the  form  of  the  disturbing  function  considered  here  is  trun¬ 
cated.  The  resonances  implied  from  equations  3.13a  through  3.13d  cause  behavior  that 
exceeds  the  limitations  of  the  second  order  theory.  The  values  u  and  fj,  are  functions  of  the 
systems  eigenvectors  and  have  the  form 


N 

(3.14a) 

N 

(3.14b) 

j=i 


Once  again,  these  relationships  imply  that  the  vector  components  and  q  consist  of 

two  parts  -  the  proper  component  and  the  forced  component: 


CpSiuWp  +  ef  sint47/ 

(3.15a) 

Bp  cos  Wp  +  Bf  cos  Wf 

(3.15b) 

Ipsinilp  +  IfsinQf 

(3.15c) 

Ip  cos  ^lp  +  If  cos  fi/. 

(3.15d) 

The  solutions  for  h,  k,  p,  and  q  can  be  used  to  determine  the  components  of  the  os¬ 
culating  orbital  elementi^^  the  orbital  elements  associated  with  the  observable  dynamics  of 
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the  particle.  The  relationship  between  these  components  are  best  thought  about  as  vector 
components.  This  relationship  can  be  depicted  graphically  through  the  use  of  a  vector 
diagram  such  as  that  shown  in  figure  3-1.  This  particular  diagram  illustrates  the  relation¬ 
ship  between  the  osculating,  forced,  and  proper  elements  of  eccentricity  and  longitude  of 
pericenter.  A  similar  diagram  can  be  constructed  for  inclination  and  longitude  of  ascending 
node  as  well. 

To  understand  the  vector  illustration,  two  important  points  must  be  understood.  First, 
the  forced  elements  are  a  function  of  the  perturbing  planets’  masses  and  osculating  elements. 
Since  the  osculating  elements  evolve,  the  forced  elements  are  also  a  function  of  time.  How¬ 
ever,  for  a  given  planetary  system  at  a  given  time,  the  forced  elements  are  a  function  of 
semi-major  axis  alone,  and  therefore  the  vector  diagram  depicts  the  relationship  between 
components  of  the  orbital  elements  for  a  particular  value  of  semi-major  axis.  Second, 
though  the  forced  elements  are  constant  for  the  particular  value  of  semi-major  axis,  all  of 
the  proper  elements  are  not  that  simple.  Though  the  proper  inclination  and  eccentricity 
may  be  constant  for  a  particular  semi-major  axis  value,  the  proper  pericenter  and  node  are 
actually  a  distribution  of  values  over  27r  resulting  from  the  randomization  of  these  elements 
after  initial  creation. 

The  forced  orbital  elements  represent  the  influence  of  the  perturbation  by  the  sec¬ 
ondary  masses  and  are  the  result  of  performing  the  calculations  outlined  in  this  section. 
By  determining  the  values  for  the  forced  orbital  elements,  the  instantaneous  trajectory  of 
a  particle  can  be  determined  provided  the  proper  elements  can  be  determined  from  the 
initial  conditions.  Prom  the  relationships  shown  in  equations  3.15a  through  3.15d,  and  the 
solutions  shown  in  equations  3.13a  through  3.13d,  the  forced  elements  can  be  expressed  as 


e/ 

=  \Jhl  +  kl 

(3.16a) 

If 

=  \/^o  +  9o 

(3.16b) 

Wf 

=  arctan  (ho/ko) 

(3.16c) 

=  arctan  (po/go) 

(3.16d) 
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Figure  3-2:  Forced  orbital  elements  as  a  function  of  semi-major  axis  for  a  massless  test 
particle  perturbed  by  the  foiu-  Jovian  planets.  The  filled  circles  represent  the  osculating 
elements  of  the  planets,  a)  forced  eccentricity  b)  forced  longitude  of  pericenter  c)  forced 
inclination  d)  forced  longitude  of  ascending  node. 

Examples  of  the  calculated  forced  orbital  elements  for  a  massless  test  particle  are  shown  in 
figure  3-2,  for  the  J2000  orbital  elements  given  in  the  appendix.  This  set  of  forced  elements 
is  based  on  perturbations  of  the  four  Jovian  planets  in  the  solar  system.  The  osculating 
elements  associated  with  each  planet  are  shown  as  filled  circles  on  the  plots.  An  important 
result  of  the  calculations  is  that  the  forced  parameters  are  equal  to  the  osculating  elements 
of  the  perturbing  planets  at  the  semi-major  axis  of  the  planets.  Further  discussion  of  the 
physical  implications  of  the  results  from  secular  perturbations  is  presented  in  the  following 
chapter.  The  discussion  will  emphasize  the  role  the  forced  inclination  plays  on  dust  particle 
orbits  within  a  disk  and  the  resulting  overall  morphology. 

Though  orbital  parameters  define  a  particle’s  trajectory,  a  trajectory  defined  by  the 
forced  orbital  parameters  alone  does  not  define  a  true  physical  particle  orbit.  As  stated 
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above,  the  forced  elements  are  only  one  part  of  the  story.  However,  by  examining  the 
forced  orbital  elements  and  trajectories  determined  by  the  forced  elements  alone,  insight 
to  a  system’s  dynamics  can  be  gleaned.  If  the  forced  elements  are  defined  as  a  function 
of  semi-major  axis  and  trajectories  based  on  the  forced  elements  are  plotted  for  a  range  of 
semi-major  axes,  the  planes  about  which  the  particles  will  precess  can  be  visualized. 

The  following  chapter  discusses  the  role  played  by  the  calculated  proper  and  forced 
orbital  elements  on  the  morphology  of  a  simulated  disk. 


CHAPTER  4 

RESULTS  FROM  PERTURBATION  THEORY 


The  previous  chapter  discussed  the  formalism  of  secular  perturbation  theory.  Within 
the  solutions  presented  in  equations  3.15a  through  3.15d  the  behavior  of  the  osculating 
orbital  elements  can  be  gleaned.  However,  the  goal  of  this  research  is  to  understand,  both 
quantitatively  and  qualitatively,  the  role  that  the  perturbations,  and  hence  the  planets,  play 
in  defining  the  morphology  of  a  circumstellar  disk.  The  goal  of  this  chapter  is  to  develop  a 
more  intuitive  understanding  about  how  the  secular  perturbations  affect  a  disk,  specifically 
observed  warping.  This  is  pursued  by  first  studying  the  details  of  the  results  presented 
in  the  previous  chapter  and  then  by  considering  the  results  in  the  context  of  the  angular 
momentum  of  dust  particle  and  planetary  orbits. 

4.1  Effects  of  Forced  Elements  on  Disk  Morphology 
As  shown  in  the  previous  chapter,  the  forced  orbital  elements  for  a  test  particle  are 
determined  by  solving  the  difierential  equations  describing  the  long  term  gravitational  per¬ 
turbations  from  massive  bodies.  It  was  also  shown  that  the  forced  elements  are  a  function 
of  semi-major  axis.  This  dependency  on  semi-major  axis  results  in  several  important  char¬ 
acteristics  that  influence  the  morphology  of  a  model  disk  that  is  based  on  these  parameters. 
An  example  of  the  forced  orbital  parameters  for  a  massless  test  particle  under  the  influence 
of  two  planets  is  illustrated  in  figure  4-1.  This  particular  calculation  considered  perturba¬ 
tions  by  a  planetary  configuration  based  on  Jupiter  and  Saturn.  Figure  4-1  illustrates  the 
dependence  of  the  forced  elements  as  a  function  of  semi-major  axis.  Of  particular  inter¬ 
est  for  this  research  are  the  characteristics  of  the  forced  elements  that  cause  the  warping 
of  a  circumstellar  disk.  This  chapter  will  discuss  in  physical  terms  how  the  variations  in 
the  forced  parameters,  particularly  the  forced  inclination,  are  manifested  into  a  warped  disk. 
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Figure  4-1:  Forced  orbital  eleraents  as  a  function  of  semi-major  axis  for  massless  test 
particles  perturbed  by  a  two  planet  system.  The  filled  circles  represent  the  osculating 
orbital  elements  of  the  perturbing  planets. 
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4.2  Solutions  from  Secular  Perturbation  Theory 

The  solutions  to  Lagrange’s  equations  of  planetary  motion  in  secular  perturbation 
theory  show  that  a  particle’s  osculating  orbital  elements  can  be  broken  down  into  two 
components:  the  forced  and  proper  orbital  parameters.  The  forced  elements  are  based  on 
the  pertiurbations  by  the  other  masses  within  the  system  and  the  proper  orbital  parameters 
are  based  on  the  problem’s  boundary  conditions.  More  specifically,  the  proper  elements 
are  related  to  the  orbital  elemenents  of  the  massless  test  particles  at  the  moment  of  their 
creation.  In  general,  this  is  determined  by  considering  the  origin  of  the  particles.  If  the  test 
particles  are  assumed  to  come  from  the  break  up  of  larger  bodies,  the  proper  elements  are 
related  to  the  characteristics  of  those  original  bodies. 

In  evolved  circumstellar  disks,  the  type  considered  in  this  research,  the  material  is 
not  primordial  but  is  indeed  created  through  the  break-up  of  leurger  parent  bodies,  either 
asteroids  or  comets.  Evolved  disks  are  disk  systems  in  which  planetesimal  formation  has 
already  occurred  or  is  in  the  process  of  occurring.  Furthermore,  radiation  effects  have 
already  had  the  opportunity  to  remove  material  via  PR  drag  or  radiation  pressure.  In 
these  types  of  systems,  the  characteristics  of  the  parent  bodies’  orbits  determine  the  proper 
orbital  elements  for  the  dust  particles.  This  is  based  on  the  idea  that  the  most  probable 
collision  in  which  a  large  body  is  broken  up,  is  one  in  which  the  colliding  particle  has  just 
enough  mass  and  energy  to  do  so.  Therefore,  the  relative  velocity  between  the  resulting 
fragments  is  low.  This  results  in  the  fragments  and  dust  having  essentially  the  same  orbital 
parameters  as  the  original  parent  body  for  inclination,  eccentricity,  and  semi-major  axis. 
However,  since  the  particles  have  a  wide  range  of  sizes  they  have  differing  values  of  jS.  This 
essentially  causes  each  size  of  particle  to  exist  in  a  unique  potential.  In  other  words,  it 
makes  each  particle  see  a  star  of  a  different  mass.  The  net  effect  of  these  differences  is 
that  the  particles  have  slightly  different  orbits  and  in  turn  disperse.  The  way  in  which  this 
is  incorporated  into  disk  models  is  by  randomizing  the  proper  pericenter  and  longitude  of 
ascending  node  over  27r  radians. 

In  the  models  created  in  this  research,  the  proper  orbital  elements  of  the  source  particles 
are  based  on  distributions  consistent  with  the  main  belt  asteroids  within  the  solar  system. 
Though  it  is  conceivable  that  multiple  reservoirs  of  source  bodies  could  exists  at  different 
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Figure  4-2:  Proper  precessional  rate  as  a  function  of  semi-major  axis.  The  horizontal 
lines  indicate  the  eigenfrequencies  associated  with  the  system:  the  solid  lines  are  eccen¬ 
tricity /pericenter  eigenfrequencies  and  the  dashed  axe  inclination/node  eigenfrequencies. 
The  intersections  between  the  eigenfrequencies  and  precessional  frequency  determine  the 
locations  of  secular  resonances. 

radii  from  a  star  and  have  distributions  of  orbital  elements  that  are  significantly  different, 
this  alone  would  not  produce  warping  in  the  disk.  It  will  be  shown  in  this  chapter  that 
varying  the  distributions  of  proper  inclination  affects  the  morphology  of  a  disk  by  changing 
the  observed  opening  angle  of  the  disk,  but  not  by  creating  a  warp. 

On  the  other  hand,  the  variation  in  forced  inclination  with  semi-major  axis  does  cause 
warping  of  the  disk.  If  it  is  accepted  that  the  forced  inclination  defines,  in  part,  the  plane 
about  which  particle  orbits  will  process,  this  is  obvious. 

By  examining  the  plot  of  the  forced  inclination  in  figure  4-1,  several  key  features  should 
be  noted.  The  first  are  the  regions  associated  with  the  singularity  in  the  secular  solutions; 
these  are  the  large  spikes  observed  in  the  function.  These  singularities  occur  at  locations 
where  the  system  eigenfrequencies  are  equal  to  the  local  precessional  rates  of  the  proper 
elements  and  represent  the  location  of  secular  resonances.  Recall  that  eigenfrequencies  are 
characteristics  of  a  dynamic  system  and  define  the  frequencies  at  which  secular  resonances 
occur.  An  example  of  how  this  is  determined  is  shown  in  figure  4-2.  This  plot  shows  the 
proper  precessional  rate.  A,  as  a  function  of  semi-major  axis.  Plotted  as  horizontal  lines  are 
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Figure  4-3;  Forced  inclination  as  a  function  of  semi-major  axis  for  massless  test  particles 
perturbed  by  a  two,  three,  and  four  planet  system.  The  numbers  indicate  the  inclination 
as  a  function  of  semi-major  axis  for  the  perturbing  planets. 

the  the  system’s  eigenfrequencies.  As  seen  in  the  equations  3.13a-3.13d,  singularities  exist 
where  the  eigenfrequencies  intersect  with  the  proper  precessional  rate.  The  second  impor¬ 
tant  characteristic  is  that  the  values  of  the  forced  elements  at  the  semi-major  axis  values 
associated  with  the  perturbing  planets  are  equal  to  the  corresponding  planets’  osculating 
orbital  elements.  This  is  illustrated  in  the  lower  left  panel  of  figure  4-1  by  plotting  the 
inclinations  of  the  perturbing  planets  as  a  function  of  their  respective  semi-major  axes.  As 
will  be  discussed  in  the  next  section,  this  correlation  between  the  forced  parameters  and  the 
parameters  associated  with  the  perturbing  planets  seems  reasonable  based  on  arguments 
made  in  terms  of  angular  momentum. 

The  same  characteristics  appear  when  the  analysis  is  performed  for  three  or  more 
planets.  Plots  of  the  forced  inclination  for  two,  three,  and  four  planets  are  shown  in  figure 
4-3.  The  number  of  resonances  increase  as  the  number  of  planets  increase  because  of  the 
additional  eigenfrequencies.  Important  for  the  following  section  is  the  fact  that  the  forced 
parameters  are  still  equal  to  the  osculating  orbital  parameters  of  the  planets  at  the  semi¬ 
major  axis  values  corresponding  to  the  semi-major  axis  of  the  planets. 
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4.3  Defining  a  Section  of  Disk 
4.3.1  General  Characteristics 

The  solutions  for  h,  k,  p,  and  q  can  be  used  to  determine  the  components  of  the 
osculating  orbital  elements  -  the  orbital  elements  associated  with  the  observable  dynamics 
of  the  particle.  The  relationship  between  these  components  is  best  considered  as  vector 
components.  The  relationship  between  components  can  be  depicted  graphically  through 
the  use  of  a  vector  diagram  such  as  that  shown  in  figure  4-4.  This  particular  diagram 
illustrates  the  relationship  between  the  osculating,  forced,  and  proper  elements  of  inclination 
and  longitude  of  ascending  node.  This  figure  is  similar  to  the  diagram  constructed  for 
eccentricity  and  pericenter  as  shown  in  the  previous  chapter. 

To  understand  the  vector  illustration,  two  important  points  must  be  understood.  First, 
for  a  given  planetary  system  at  a  given  time  the  forced  elements  are  a  function  of  semi-major 
axis  alone,  and  therefore  the  vector  diagram  depicts  the  relationship  between  components  of 
the  orbital  elements  for  that  particular  value  of  semi-major  axis.  Second,  though  the  forced 
elements  are  single  valued  for  a  particular  value  of  semi-major  axis,  the  proper  elements  are 
not.  The  proper  pericenter  and  node,  for  a  particular  semi-major  axis  value,  are  actually 
a  distribution  of  values  over  27r  resulting  from  the  randomization  of  these  elements  after 
initial  creation.  As  mentioned  earlier,  this  results  from  phase  mixing  -  the  dispersal  of 
particles  within  an  ensemble  due  to  slight  differences  in  orbital  periodicity. 

Figure  4-5  shows  the  relationship  between  the  components  of  the  forced  and  proper 
orbital  elements  and  the  section  of  disk  which  they  define.  Because  the  proper  pericenter 
and  node  are  randomized  over  27r,  particles  at  a  given  value  of  semi-major  axis  will  create  a 
distribution  of  orbits,  each  with  slightly  different  osculating  elements.  This  distribution  will 
result  in  a  torus  that  is  responsible  for  defining  that  particular  section  of  disk.  These  plots 
show  the  distribution  of  orbits  having  the  same  semi-major  axis,  and  thus  the  same  forced 
eccentricity  and  the  same  forced  pericenter.  The  top-left  panel  has  a  forced  inclination  of 
zero  while  the  other  panels  show  disks  with  a  nonzero  forced  inclination.  The  filled  circle 
represents  the  geometric  center  of  the  torus  and  and  the  empty  circle  represents  the  central 


star. 


38 


X-Y  plone  X-Y  plane 


-1,5 

; . 

1  ....  t  ....  1  ....  1  ....  n 

-1.5 

~r  ■  .  .  ■  1  ■  ■  ■  .  1  ■  .  ■  ■  1  ■  >  ■  ■  t  ■  .  ■  .  »  ■  .  ■  i-n 

1.5  -1.0 

-0.5  0.0  0.5  1.0  1.5 

X~-Z  plane 

1.5  -1.0  -0.5  0.0  0.5  1.0  1.5 

Y-Z  plane 

1.5 

j  .  1  '  '  1  ■  ' 

:  (c) 

'  '  I  '  '  M  '  *  '  '  1  ■  ■  *  ■  1  '  '  '  '  L 

1.5 

J  ’  '  '  '  1  '  '  '  '  I  •  ■  ■  ‘  I  -  .  .  ■  1  ■  ■  ■  ■  1  •  .  /-r-|_ 

■  (d)  ; 

1.0 

- 

1.0 

“ 

0.5 

- 

—  ; 

0.5 

0.0 

0.0 

-  •.v'--"''-  ...  if- 

-0.5 

-0.5 

-1.0 

-1.0 

-  ^ 

-1.5 

■  i  ■  ■ 

..  1  ....  . . .  .  .  .  .  i~ 

-1.5 

h 

. . .  . . . . .  r 

-1.5  -1.0  -0.5  0.0  0.5  1.0  1.5  -1.5  -1.0  -0.5  0.0  0.5  1.0  1.5 


Figure  4-5:  The  distribution  of  test  particle  orbital  trajectories  representing  the  influence 
of  forced  elements.  Panel  a  shows  the  torus  resulting  from  a  disk  with  a  nonzero  forced  and 
proper  eccentricity.  The  Other  panels  show  different  views  of  (a)  the  same  disk  but  with 
nonzero  values  of  forced  and  proper  inclination  with  respect  to  the  X-Y  plane.  The  forced 
eccentricity  determines  the  ofiset  of  the  disks  center  of  symmetry  from  the  central  star  and 
the  proper  eccentricity  determines  the  width  of  the  tauri.  Similarly  the  forced  inclination 
determines  the  angle  of  the  orbital  plane  to  a  reference  plane  and  the  proper  inclination 
determines  the  disk  opening  angle. 
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Figure  4-6:  Visualization  of  a  section  of  disk.  The  offset  of  the  disk  section’s  center  of 
symmetry  is  ~  ae/.  The  angular  thickness  is  ^  2ip  and  the  angle  of  the  section’s  symmetry 
plane  with  respect  to  the  system’s  reference  plane  is  ~  iforced- 

All  of  the  characteristics  of  the  torus  representing  the  distribution  of  dust  particle  orbits 
can  be  related  to  the  forced  and  proper  elements.  Generally  speaking,  the  proper  elements 
define  the  geometric  size  of  the  torus  while  the  forced  elements  define  the  orientation  of 
the  torus  with  respect  to  the  central  star  and  an  associated  reference  plane  and  direction  - 
such  as  the  ecliptic  plane  and  the  vernal  equinox  for  our  solar  system. 

The  dimensions  of  the  torus  are  related  to  the  proper  elements  of  eccentricity  and 
inclination.  The  inner  radius  of  the  torus  is  equal  to  a(l  —  Sp)  and  the  outer  radius  is  equal 
to  a{l  +  ep).  The  center  of  symmetry  is  offset  from  the  central  star  by  a  distance  aej  in  the 
direction  opposite  to  Wf  (Dermott  et  al.  1985;  Dermott  et  al.  1998). 

Likewise,  the  proper  inclination  determines  the  angular  thickness  of  the  disk.  As  has 
been  shown  for  the  solar  system  dust  bands  (Dermott  et  al.  2001),  the  angular  thickness 
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Figure  4-7:  Proper  precessional  rates  and  forced  inclination  for  a  system  modeled  after 
Jupiter  and  Saturn.  The  top  illustrations  show  how  the  precessional  rate/  inclination  vary 
with  radial  distance  from  the  central  star  and  with  an  increase  in  the  mass  and  semi-major 
axes  of  the  planets.  The  bottom  plots  show  the  same  for  the  first  and  last  values  of  mass 
and  semi-major  axes  scale  factor. 


~  2Ip.  The  forced  inclination  determines  the  plane  of  symmetry  or  the  disk  section’s  mid¬ 
plane.  This  is  illustrated  in  panel  (d)  of  figure  4-5  which  shows  a  cross-sectional  view  of 
the  torus  representing  the  contributions  of  dust  particles  with  a  common  semi-major  axis 
to  the  overall  disk.  An  actual  disk  could  be  though  of  as  a  superposition  of  many  such 
tori,  each  characterized  by  a  set  of  forced  parameters  associated  with  a  particular  value 
of  semi-major  axis.  Figure  4-6  shows  another  view  of  a  disk  section  and  illustrates  the 
relationships  between  the  geometry  and  the  forced  and  proper  elements. 


4.3.2  Determining  the  Radial  Extent  of  the  Influence  of  Perturbing  Planets 

Analysis  of  the  forced  orbital  elements  determined  from  secular  perturbation  theory 
show  an  invariant  value  is  approached  at  large  vales  of  semi-major  axis.  The  region  through 
which  the  elements  begin  to  approach  this  invariant  value  define  the  extent  of  observable 
effects  from  the  secular  perturbations. 
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Figure  4-7  provides  an  illustration  of  how  the  range  of  influence  for  a  given  planetary 
configuration  is  altered  by  modifying  both  the  planetary  masses  and  the  semi-major  axis 
of  the  planets  by  an  increasing  scale  factor.  The  top  surface  plots  show  the  changes  to 
the  proper  precessional  rate,  A,  and  the  forced  inclination,  If.  The  diagram  shows  the 
amplitude  of  A  and  If  with  semi-major  axis  and  increasing  values  of  the  mass  and  semi¬ 
major  axis  scale  factor.  The  bottom  plots  show  a  two  dimensional  version  of  the  same,  with 
only  the  first  and  last  values  of  the  scale  factor  plotted. 

The  increase  in  mass  and  the  increase  in  semi-major  axis  have  unique  effects  on  the 
elements.  The  semi-major  axis  scale  factor  determines  the  location  of  the  important  features 
in  the  plots  including  the  locations  of  the  singularities  as  well  as  the  location  of  the  invariant 
plane  for  the  forced  elements. 

Increasing  the  mass  scale  factor  has  a  different  effect  in  that  it  extends  the  region  over 
which  an  individual  planet  exerts  an  influence.  For  example,  the  width  of  the  singularities 
associated  with  the  proper  precessional  rates  are  widened  as  the  masses  of  the  planets 
increase.  Since  the  proper  precessional  rate  also  determines  the  timescales  of  secular  effects, 
this  effect  highlights  how  more  massive  systems  have  shorter  secular  timescales. 

An  important  characteristic  of  the  forced  elements  is  that  they  match  the  osculating 
elements  of  the  perturbing  planets  when  the  particle’s  semi-major  axis  matches  that  of  the 
corresponding  planet.  As  can  be  seen  in  the  plot  of  inchnation,  a  small  plateau  exists  in  the 
vicinity  of  the  planet.  For  more  massive  planets,  this  plateau  is  increased  in  its  semi-major 
axis  range.  This  implies  that  a  more  massive  planet  will  impose  its  orbital  elements  on  a 
test  particle  over  a  larger  region  of  a  disk. 

By  developing  intuition  on  how  these  elements  behave  when  scaling  parameters  such  as 
semi-major  axis  and  planetary  mass,  and  by  understanding  how  these  elements  ultimately 
affect  the  morphology  of  a  disk,  developing  an  initial  planetary  configuration  for  a  model 
disk  becomes  much  less  of  a  process  based  solely  on  trial  and  error. 

4.4  Planes  Defined  by  the  Angular  Momentum  Vector 

As  stated  in  the  previous  section,  a  warped  disk  is  the  result  of  the  forced  inclination 
varying  as  a  function  of  semi-major  axis.  This  variation  in  forced  inclination,  as  is  the 
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Figure  4-8:  Planes  defined  by  angular  momentum  of  individual  planets  in  a  multi-planet 
system.  The  reference  plane  represents  the  invariant  plane;  that  plane  perpendicular  to  the 
system’s  total  angular  momentum  vector. 

case  with  the  other  orbital  parameters  as  well,  is  the  result  of  the  gravitational  perturba¬ 
tions  associated  with  the  presence  of  planets.  Another  way  to  consider  why  this  occurs 
is  by  considering  the  planes  associated  with  angular  momentum  vectors;  both  the  angular 
momentum  associated  with  each  planet  and  the  total  angular  momentum  of  the  system. 

Warps  in  circumstellar  disks  have  been  indicated  in  both  observations  of  real  disks  and 
theoretical  models.  In  both  cases,  the  warping  has  been  described  as  the  edge  on  view  of 
a  multi-component  disk  in  which  the  components  are  inclined  to  one  another.  The  origin 
of  these  multiple  components  are  generally  attributed  to  the  presence  of  planets.  One  of 
the  more  thorough  analyses  argues  that  the  multiple  component  configuration  results  from 
a  single  planet  having  an  orbit  inclined  to  the  plane  of  the  disk.  The  general  idea  is  that 
through  interaction  between  the  planet  and  the  disk  material,  the  dust  becomes,  over  time, 
coplanar  with  the  orbital  plane  of  the  planet.  Furthermore,  the  presence  of  the  planets 
cause  differential  rotation  throughout  the  disk.  If  the  system  is  relatively  young,  material 
at  outer  radii,  whose  orbits  process  more  slowly,  has  not  yet  become  coplanar.  This  scenario 
results  in  two  planes  of  symmetry  and  satisfies  the  premise  for  the  multi-component  theory. 

However,  there  exists  a  more  evolved  system  in  which  a  warp  is  present  as  well:  the 
solar  system  dust  bands.  Therefore,  there  must  be  another  mechanism  to  create  the  multiple 
symmetry  planes  or  components.  If  one  considers  the  angular  momentum  associated  with 
individual  planetary  orbits  as  well  as  with  the  entire  multi-planet  system,  the  multiple 
components  naturally  arise. 
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Arguments  made  considering  the  influence  of  a  single  planet  on  a  circumstellar  disk, 
describe  a  situation  in  which,  through  torquing,  particle  orbits  will  slowly  become  copla- 
nar  with  the  plane  defined  as  being  perpendicular  to  the  angular  momentum  vector  of  the 
planet’s  orbit.  This  can  be  extended  to,  but  becomes  more  complicated  for,  a  multi-planet 
system  as  each  planetary  orbit  has  associated  with  it  an  angular  momentum  vector.  Fur¬ 
thermore,  there  exists  a  total  angular  momentum  vector  for  the  system  which  defines  the 
plane  at  large  distances  from  the  system,  where  the  radii  of  the  individual  masses  are  com¬ 
paratively  small.  To  state  more  succinctly,  multiple  planes  or  components  will  naturally 
arise  through  the  angular  momentum  in  a  multi-planet  system. 

Determination  of  the  orbital  angular  momentum  is  based  on  the  mass  of  the  planet 
and  the  planet’s  orbital  parameters:  semi-major  axis,  mean  rate,  inclination,  and  longitude 
of  ascending  node.  Therefore,  two  bodies  with  angular  momentum  vectors  oriented  in 
the  same  direction  will  have,  at  the  very  least,  the  same  orbital  plane.  Using  the  orbital 
parameters  of  a  three  planet  system,  the  planes  perpendicular  to  the  individual  angular 
momentum  vectors  have  been  calculated  and  are  shown  in  figure  4-8.  The  reference  plane 
is  the  invariant  plane  associated  with  the  system’s  total  angular  momentum  vector.  The 
values  used  in  this  plot  are  the  same  as  those  used  in  creating  the  middle  panel  of  figure 
4-3. 

By  carefully  considering  the  information  in  figure  4-8  and  the  corresponding  panel  in 
figure  4-3.  relationships  can  be  discerned.  First,  if  the  forced  orbital  elements  match  those 
of  one  of  the  perturbing  planets  at  the  same  semi-major  axis  of  that  planet,  the  test  particle 
orbits  necessarily  will  be  precessing  around  a  plane  that  is  in  fact  perpendicular  to  that 
specific  planet’s  orbital  angular  momentum  vector;  i.e.  one  of  the  planes  shown  in  figure 
4-8.  Second,  on  large  scales,  the  particle  orbits  will  process  around  an  invariant  plane  that 
is  perpendicular  to  the  systems  total  angular  momentum  vector;  i.e.  the  reference  plane 
shown  in  figure  4-8. 

As  should  be  expected,  other  permutations  of  this  model  showed  that  the  number 
of  components  increases  with  the  number  of  planets.  This  is  consistent  with  the  results 
illustrated  in  figiure  4-9  which  show  how  each  planet  has  associated  with  it  a  plane  defined 
by  its  angular  momentum  vector. 
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It  should  be  noted  that,  though  these  models  have  been  designed  to  highlight  the  po¬ 
tential  for  warping,  it  could  easily  be  imagined  that  a  system  could  be  composed  of  many 
planets  and  still  show  little  sign  of  a  multiple  component  configuration  if  the  corresponding 
nodes  were  not  aligned  and  the  mutual  inclination  was  low. 

4.5  Model  Disks 

Simulations  of  model  dust  disks  were  created  to  support  the  suggestions  made  within 
the  previous  sections  that  the  variation  of  forced  inclination  with  semi-major  axis  is  capa¬ 
ble  of  producing  observable  warping  of  the  type  identified  in  the  Beta  Pictoris  disk.  These 
simulations  involved  first  calculating  the  forced  orbital  parameters  for  a  particular  plane¬ 
tary  configuration.  Second,  using  assumed  values  for  the  proper  orbital  elements,  a  model 
dust  disk  was  generated.  The  resulting  models,  shown  in  figure  4-9  are  three  dimensional 
mappings  of  surface  area.  The  primary  significance  of  these  diagrams  is  in  showing  the 
symmetry  planes  for  a  disk  perturbed  by  a  particular  configuration  of  planets. 

Different  views  of  three  model  disks  are  shown  in  figure  4-9.  The  top  row  depicts  a 
disk  pertmbed  by  two  planets;  the  second  row  depicts  a  disk  perturbed  by  three  planets; 
and  the  bottom  row  depicts  a  disk  perturbed  by  four  planets.  The  first  colunm  shows  a 
face-on  view  of  the  model  disks.  The  second  column  shows  a  nearly  edge-on  view,  ~  3 
degrees  of  tilt  to  the  line  of  sight.  The  third  column  shows  the  same  image  presented  in 
the  second  column  but  stretched  vertical  by  a  factor  of  fomr.  Furthermore,  the  images  in 
the  third  column  have  hawl  each  vertical  cut  in  the  disk  image  normalized  to  the  peak  value 
of  that  particular  cut.  The  calculated  forced  inclinations  for  each  planetary  configuration 
used  in  developing  each  disk  model  are  shown  in  figure  4-3. 

Warping  is  clearly  seen  in  each  model,  especially  in  the  stretched  and  normalized 
images.  Within  each  image,  one  can  easily  visualize  a  configuration  of  multiple  components 
which  are  offset  to  one  another.  Furthermore,  the  inclusion  of  additional  planets  seem  to 
add  more  components  and  additional  structure  to  the  disk.  As  was  explained  in  an  earlier 
section,  this  seems  intuitive  if  one  considers  this  problem  in  terms  of  angular  momentum. 

Comparisons  can  also  be  made  between  the  four  planet  simulations  and  the  plot  of 
angular  momentum  planes  illustrated  in  figure  4-8,  which  are  both  based  on  the  same 


-100  -50  0  50  100 

Figure  4-10:  Planes  defined  by  angular  momentum  of  individual  planets  in  a  multi-planet 
system.  The  reference  plane  represents  the  invariant  plane;  that  plane  perpendicular  to  the 
system’s  total  angular  momentum  vector. 
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orbital  parameters.  In  the  simulations  shown  in  figure  4-9,  multiple  components  can  be 
identified.  Furthermore,  the  number  of  components  increases  with  the  number  of  planets. 
This  is  consistent  with  the  results  illustrated  in  figure  4-10,  which  show  how  each  planet 
has  associated  with  it  a  plane  defined  by  its  angular  momentum.  This  figure  is  a  simulated 
disk  with  orbital  trajectories  of  the  perturbing  planets  superimposed  over  the  disk.  This 
illustrates  how  the  orbital  planes  are  coplanar  to  the  major  symmetry  planes  of  the  model 
dust  disk. 


4.6  Inclination  Excitation 

4.6.1  Introduction 

The  previous  sections  have  discussed  how  the  mid-plane  of  a  circumstellar  disk  can  be 
warped  due  to  the  presence  of  a  multi-planet  system.  One  of  the  requirements  for  producing 
the  mid-plane  warping  is  a  relatively  large  mutual  inclination  between  the  planets  within  the 
system.  In  the  first  chapter,  it  was  asserted  that  previous  models  of  the  disk  morphology 
were  questionable  due  to  their  unclear  explanation  of  exactly  how  a  single  planet  would 
become  highly  inclined  to  the  disk  from  which  it  formed.  If  the  model  of  the  disk  presented 
in  this  research  is  indeed  more  plausible,  then  the  origin  of  the  high  mutual  inclination 
between  the  constituent  planets  must  be  addressed. 

Fortunately,  theoretical  work  on  orbital  element  excitation  has  already  been  well  stud¬ 
ied.  Dermott,  Malhorta,  and  Murray  (1988)  illustrated  how  inclination  and  eccentricity 
can  be  increased  dramatically  by  the  temporary  capture  of  a  satellite  into  a  resonance.  In 
this  work  it  was  shown  that  for  two  satellites  moving  on  a  converging  orbit,  the  evolution 
that  occurs  during  capture  will  increase  the  orbital  elements.  This  is  based  on  tidal  inter¬ 
actions  that  work  to  transfer  angular  momentum  firom  the  inner  to  the  outer  body.  Chiang 
(2003)  discussed  how  this  can  occur  in  exosolar  planets  that  are  migrating  in  a  circumstellar 
disk.  Presented  in  this  section  is  summary  of  the  mathematics  behind  these  theories  and  a 
discussion  of  the  applicability  to  planetary  systems  on  a  scale  similar  to  that  of  Beta  Pictoris. 
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4.6.2  Excitation  of  Orbital  Elements  by  Capture  into  Resonance 

After  protoplanets  have  formed  and  accretion  has  stopped  due  to  the  clearing  of  a 
gap  in  the  disk,  a  protoplanet  continues  to  interact  with  disk  through  tidal  forces.  This 
interaction  can  lead  to  an  evolution  of  the  protoplanet’s  orbit  through  eccentricity  excitation 
and/or  orbital  migration  (Lin  and  Papaloizou  1993). 

The  evolution  of  the  orbital  elements  associated  with  planets  orbiting  a  central  star 
can  be  described  through  the  use  of  Lagrange’s  planetary  equations.  Lagrange’s  equations 
for  the  mean  rate,  eccentricity,  and  inclination  are  given  below: 
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where  a,  e,  7,  w,  fl,  and  A  are  the  semi-major  axis,  eccentricity,  inclination,  longitude  of 
pericenter,  longitude  of  ascending  node,  and  mean  longitude  respectively. 

The  disturbing  functions  for  m  and  m'  were  given  in  the  previous  chapter  and  were 
expressed  as 
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Prom  the  approach  of  Dermott  et  al.  (1988)  these  expressions  can  be  expanded  in  a  Fourrier 
series  and  expressed  as 


Ti  =  fi'S  cos  (j) 
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where  /x  =  Gm  and  fj.'  =  Gm'.  The  S^term  can  be  thought  of  as  the  ’strength’  of  a  particular 
term  in  the  expansion  and,  to  lowest  order,  can  be  given  as 

S  =  ZMel9ile'l«2lsl93l5'k4|  (4.5) 

a' 

where  f{a)  can  be  expressed  as  a  function  of  Laplace  coefficients  and  s  =  sin  The  values 
of  q  represent  the  commensurability  between  the  periods  of  the  planets’  orbits,  in  the  form 
oip  :p  +  q.  Using  these  variables,  the  term  (j)  is  given  as 

(j)  =  pX-  {p  +  q)X'  +  q\w  +  q2XX}'  +  q^^  +  q^^'  (4.6) 

where  the  value  q  is  given  as 

Q  =  (4-7) 

the  time  rate  of  change  of  the  mean  rate,  n,  can  be  considered  to  be  changing  due  to 
resonance  and  tidal  effects.  Using  equation  4.1a,  the  change  in  mean  rate  can  be  written  as 

n  =  cos  +ht 
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and  likewise 
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An  expression  for  5sin^,  can  be  found  by  taking  the  second  time  derivative  of  (f> 

0  =  pA  -  (p  +  q)X'  +  qiw  +  q2w'  +  qz^  +  q^^'  (4-10) 

Since  the  second  derivatives  of  the  orbital  elements  are  generally  small,  they  can  be  ignored. 
This  allows  (f>  to  be  expressed  as 

-3  3 
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where  the  values  of  g  and  F  are  given  by 

9  =  +  +  (4-12) 

and  F  =  pht  —  (p  +  q)n[  (4-13) 

where  F  is  determined  from  the  drag  forces  on  the  body.  An  expression  for  the  average 
value  of  5sin^  can  be  found  by  assuming  that  the  <  ^  >=  0  which  leads  to 

(Ssin^)  = 

The  average  values  for  /  can  be  expressed  using  the  equation 
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If  the  problem  is  constrained  to  low  inclinations  and  eccentricities 

(j)  = 

similarly,  the  expression  for  the  time  rate  of  change  of  orbital  eccentricity  can  be  written  as 


(4.17) 


(i) = -^^4 

It  needs  to  be  pointed  out  that  though  equations  4.17  and  4.18  may  seem  to  indicate 
that  higher  order  resonances  play  bigger  roles,  through  the  qi  terms,  this  is  not  the  whole 
story.  The  probability  of  capture  into  a  resonance  must  be  considered  and  this  is  smaller 
for  the  higher  order  resonances  (Dermott  et  ah,  1988). 

Equations  4.17  and  4.18  show  how  the  inclination  and  eccentricity  can  grow  if  there 
exists  tidal  interactions  and  a  commensurability  between  the  periods  of  the  orbits.  It  is 
generally  believed  that  planets  may  migrate  during  the  planet  forming  epoch  and  are  thus 
likely  to  pass  through  resonances.  This  leads  to  the  conclusion  that  planets  should,  pro¬ 
vided  they  are  able  to  interact  tidally,  experience  excitation  to  some  extent.  This  supports 
the  argument  that  a  non-mutual  inclination  should  be  expected  in  a  multi-planet  system. 
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4.7  Conclusions 

This  is  an  investigation  motivated  by  the  observations  of  mid-plane  warping  in  real 
circumstellar  disks.  Based  on  the  idea  that  an  observed  warp  is  actually  the  edge-on  pro¬ 
jection  of  a  multi  component  disk,  in  which  the  components  are  inclined  to  one  another,  the 
goal  was  to  determine  what  mechanism  was  responsible  for  creating  the  multi-component 
configuration.  It  is  suggested  here  that  multiple  planes  of  symmetry,  or  components,  nat¬ 
urally  arise  in  a  multi-planet  system  due  to  the  differences  in  angular  momentum  between 
the  constituent  planets.  Material  within  the  disk  at  the  same  radii  as  a  planet,  behaves 
in  such  a  way  as  to  process  around  the  plane  defined  by  that  planets’  angular  momentum; 
between  the  planets  the  behavior  is  more  complicated  as  the  dynamics  become  influenced 
by  a  combination  of  gravitational  effects  from  different  planets. 

This  qualitative  argument  is  complimentary  to  the  conclusions  reached  through  per¬ 
turbation  theory.  The  results  from  secular  perturbation  theory  illustrate  how  the  forced 
orbital  elements  of  the  test  particles  are  equal  to  the  osculating  elements  of  the  perturbing 
planet  when  the  orbits  have  a  common  semi-major  axis.  Additionally,  the  mass  of  the 
planets  play  a  role  in  defining  the  region  over  which  they  impose  their  orbital  parameters 
on  a  near-by  particle. 


CHAPTER  5 

PRIOR  OBSERVATIONS  AND  MODELING  OF  BETA  PICTORIS 


Though  the  observation  of  individual  planets  is  still  out  of  reach,  recent  observations 
of  circumstellar  disks  have  provided  another  possibility  for  planetary  detection.  The  most 
well  studied  of  the  circumstellar  disks  is  that  of  Beta  Pictoris.  Beta  Pictoris  is  a  ~  20 
Myr  old  main  sequence  star  with  a  disk  that  is  large  (extending  out  to  ~1000  AU  from 
the  central  star)  ,  bright,  and  viewed  nearly  edge-on.  Even  more  interesting  than  the  size 
and  brightness  of  the  Beta  Pictoris  disk  are  the  asymmetries  and  features  observed  in  the 
morphology  of  the  disk  -  the  most  interesting  being  the  warping  of  the  inner  disk. 

This  chapter  will  provide  a  general  review  of  the  observations  made  of  scattered  light 
and  thermal  emissions  originating  from  the  Beta  Pictoris  circumstellar  disk.  The  major 
asymmetries  observed  in  the  disk  on  both  large  and  small  scales  are  discussed  and  quanti¬ 
fied  when  possible.  The  stellar  properties  of  Beta  Pictoris,  as  determined  through  analysis 
of  Hipparcos  data,  are  also  reviewed.  These  properties  are  important  in  determining  char¬ 
acteristics  like  age  and  luminosity  which  play  an  important  role  in  the  dynamics  of  the 
material  within  the  disk. 

Past  approaches  to  modeling  the  Beta  Pictoris  circumstellar  disk  are  also  discussed. 
The  two  primary  approaches  of  modeling,  a  close  stellar  encounter  or  the  presence  of  a 
single  planet  on  an  inclined  orbit,  are  reviewed  and  their  limitations  are  discussed. 

The  goal  of  this  chapter  is  to  provide  the  observational  constraints  that  will  be  applied 
to  the  simulated  circumstellar  disks  to  be  discussed  in  the  following  chapters. 

5.1  The  Discovery  of  the  Beta  Pictoris  Disk 

In  1984  Aumann  presented  a  sample  of  stars  with  infrared  fluxes  measured  at  12 
and  25  fira.  from  the  Infrared  Astronomical  Survey  (IRAS).  Approximately  18  percent  of 
the  280  stars  in  the  survey  showed  fluxes  in  the  infrared  beyond  what  was  appropriate 
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Figure  5-1:  60  nm  point  source  image  of  Beta  Pictoris  acquired  during  the  IRAS  survey 
(Aumann,  1984) 

for  an  isothermal  blackbody  (Aumann,  1984).  The  infrared  excess  was  interpreted  as  the 
superposition  of  the  blackbody  curves  of  the  star  and  the  dust  component  or  components 
each  with  a  different  effective  temperature.  Aumann  further  suggested  that  the  infrared 
excess  was  the  result  of  the  star  being  siirrounded  by  a  cloud  of  cold,  solid  material  possibly 
in  the  form  of  a  circinnstellar  disk  (Aumann,  1984b). 

The  stars  possessing  the  brightest  infrared  excess  were  Vega  (a-Lyrae),  Fomalhaut, 
Epsilon  Eridani,  and  Beta  Pictoris.  These  stars  are  at  a  distance  of  7.75,  7.69,  3.22,  and 
19.28  pc  respectively.  Of  the  three.  Beta  Pictoris,  despite  being  the  farthest,  possessed  the 
brightest  excess  (Aumann,  1985).  Based  on  the  assertion  that  a  significant  fraction  of  stars 
(approximately  30%  of  all  F  and  G  type  stars)  in  the  solar  neighborhood  were  surrounded 
by  circumstellar  disks  and  the  fact  that  Beta  Pictoris  was  the  best  example,  Beta  Pictoris 
became  a  very  well  studied  object.  The  60  fj.m  image  of  Beta  Pictoris  taken  during  the 
IRAS  survey  is  shown  in  figure  5-1. 

When  the  disk  was  first  imaged  in  scattered  light  by  Smith  and  Terrile  (1984),  the 
observation  showed  a  nearly  edge  on  disk  that  was  resolved  out  to  25  arcsec  (~  500  AU). 
Later,  higher  resolution  observations  of  the  scattered  light  and  the  thermal  emission  from 
the  dust  disk  showed  structures  and  asymmetries  that  were  argued  to  be  related  to  the 
presence  of  a  planet  or  a  planetary  system. 

The  remainder  of  this  chapter  will  review  many  of  the  observations  made  describing  the 
featmes,  structures,  and  asymmetries  of  the  Beta  Pictoris  dust  disk.  The  following  section 
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reviews  the  stellar  properties  of  Beta  Pictoris  determined  from  Hipparcos  measurements 
and  analysis  of  the  Beta  Pictoris  moving  group.  The  next  few  sections  will  discuss  the 
specific  attributes  of  the  disk  and  how  they  are  determined  through  different  observational 
methods.  Certain  features  are  seen  only  at  certain  wavelengths  due  to  the  difierences  in 
the  sources  of  light  -  i.e.  scattering  of  optical  light  verses  the  thermal  emission  of  infirared 
radiation.  Section  5.3  reviews  and  highlights  first  the  scattered  light  images  and  then  the 
thermal  emission  observations.  Also  discussed  are  spectroscopic  observations  and  measure¬ 
ments  of  the  photometric  variations  of  Beta  Pictoris  and  the  associated  implications  to  the 
structure  of  the  dust  disk. 


5.2  Stellar  Properties 

Based  on  data  gained  from  the  Hipparcos  mission,  Beta  Pictoris  was  determined  to 
be  at  a  distance  of  19.28  ±  0.19  pc  (Crifo  et  al,  1997)  .  Prom  the  distance  measurements, 
the  ratio  of  stellar  to  solar  luminosity,  L/Lq,  was  estimated  to  be  8.7.  Based  on  analysis 
of  the  H7  line  profile,  the  effective  temperature  of  Beta  Pictoris  was  estimated  to  be  8200 
K  (Crifo  et  al,  1997).  Furthermore,  using  models  of  main  sequence  and  pre-main  sequence 
tracks,  Crifo  et  al.  used  the  effective  temperature  and  the  liuninosity  to  suggest  that  Beta 
Pictoris  fell  very  close  to  the  zero  age  main  sequence  (ZAMS)  and  had  a  mass  of  1.7-1. 9 
Mq.  All  of  the  conclusions  reached  by  Crifo  et  al.  based  on  Hipparcos  data  were  consistent 
with  a  prior  spectrophotometric  study  (Alekseeva  et  al.,  1996)  of  Beta  Pictoris  which  found 
Te//=7950  K,  L  ss  8  Lq,  and  a  minimum  age  of  ~  20  Myr. 

The  accurate  determination  of  distance  determined  by  Hipparcos  changed  prior  con¬ 
clusions  based  on  an  assmned  distance  of  ~16  pc  for  Beta  Pictoris.  For  example,  a  break 
in  the  power  law  describing  the  radial  surface  brightness  distribution  determined  by  the 
observations  of  Golimowski  et  al.  (1993)  to  be  at  100  AU  for  a  distance  of  16.4  pc  became 
120  AU  for  19.28  pc.  In  addition,  the  inner  radius  of  the  central  hole  of  the  disk,  usually  es¬ 
timated  as  being  at  30  AU,  became  35  AU.  Within  this  chapter,  all  distance  measurements 
are  based  on  the  currently  accepted  distance  of  19.28  pc  as  determined  from  Hipparcos. 

A  separate  determination  of  the  age  of  Beta  Pictoris  was  made  based  on  a  study  of  the 
Beta  Pictoris  moving  group  (Navascues  et  al.,  1999).  This  study  was  based  on  determining 
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the  stellar  properties  of  several  observed  stars  deemed  dynamically  related  to  Beta  Pictoris. 
The  mean  age  of  the  group  was  estimated  as  20±10  Myr.  Assuming  aU  of  the  stars  within 
the  group  to  be  coeval,  this  age  is  applicable  to  Beta  Pictoris  as  well. 

5.3  Observations 

Following  the  initial  measurement  of  the  infrared  excess  associated  with  Beta  Pictoris 
and  the  interpretation  that  the  excess  was  associated  with  a  disk  of  material  surrounding 
the  star,  a  direct  image  of  the  disk  became  a  goal  of  the  astronomical  community.  The 
first  images  were  obtained  through  observations  of  scattered  light  using  the  technique  of 
coronagraphy.  This  approach  observed  the  dust  disk  by  blocking  out  the  light  from  the  star 
that  would,  in  normal  circumstances,  overwhelm  the  light  scattered  from  material  within 
the  disk.  As  will  be  discussed  in  the  following  sections,  early  observations  using  this  method 
showed  a  debris  disk  extending  out  nearly  48  arcsec  (~  900  AU).  Also  seen  in  these  images 
were  asymmetries  in  the  disk  that  began  the  debate  over  the  possibility  of  a  planetary 
perturbation  being  the  source  of  the  observed  structure. 

Observations  of  the  thermal  emissions  of  the  debris  disk,  particularly  around  10  ^m, 
precluded  the  need  for  techniques  such  as  coronagraphy,  as  the  intensity  of  the  thermal 
emissions  of  the  dust  was  comparable  to  that  of  the  photosphere  of  the  star.  Thermal 
images  of  the  debris  disk  were  acquired  later  and  also  showed  asymmetries.  However,  the 
features  associated  with  the  thermal  emissions  were  not  the  same  as  those  seen  in  scattered 
light.  The  asymmetries  occurred  on  different  scales  and  were,  in  some  instances,  reversed 
from  the  scattered  light  images. 

Other  observations  were  also  taken  that  provided  more  clues  to  the  composition  and 
dynamics  of  the  dust  disk.  For  instance,  spectroscopic  variations  were  used  to  mgue  the 
presence  of  cometaiy  bodies  which  were  evaporating.  Also,  measmements  of  photometric 
variations  were  used  to  validate  the  idea  that  the  the  Beta  Pictoris  disk  possessed  significant 
asymmetries  in  the  azimuthal  distribution  of  dust. 

This  section  gives  a  brief  discussion  of  several  of  the  significant  observations  made  with 
the  above  techniques  and  discusses  their  significance  in  terms  of  providing  constraints  to  a 


Figure  5-2:  890  nin  coronagraphic  image  of  Beta  Pictoris  in  scattered  light  (Smith  and 
Terrile,  1984). 

complete  dynamic  model  of  the  dust  disk. 

5.3.1  Scattered  Light 

In  1984  Smith  and  Terrile  produced  the  first  direct  image  of  the  Beta  Pictoris  disk 
using  a  coronagraph  at  890  nm.  The  image  showed  what  was  interpreted  as  an  edge-on 
disk  extending  to  25  arcsec  (~  500  AU).  Their  coronagraphic  image  is  shown  in  figure  5-2. 
The  mask  within  the  coronagraph  covered  the  central  7  arcsec  (~  135  AU)  of  the  image 
and  therefore  the  central  region  of  the  disk  was  not  imaged.  However,  it  was  inferred  from 
extrapolation  of  the  optical  thickness  that  there  exists  a  central  hole  within  the  inner  35 
AU.  The  optical  thickness  was  estimated  using  the  observed  brightness  distribution  and 
assumptions  about  the  physical  properties  and  size  distribution  of  the  scattering  particles. 
The  argument  for  the  hole  was  based  on  the  fact  that  if  the  optical  depth  continued  to 
increase  at  the  estimated  rate,  a  disk  approaching  closer  than  35  AU  would  be  inconsistent 
with  observed  attenuation  of  the  star’s  radiation.  Smith  and  Terrile  also  claimed  that  the 
disk  was  wedge-shaped,  having  a  width  that  increased  with  radius;  measurements  of  the 
angle  associated  with  the  wedge  indicated  that  the  particles  within  the  disk  were  inclined 
at  5°  or  less.  It  was  also  inferred  that  because  no  reddening  was  observed  in  multicolored 
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photometry,  the  size  of  the  emitting  particles  must  be  much  larger  that  the  wavelength  of 
the  light  observed;  A  d  ~  few  /xm.  Smith  and  Terrile  argued  that  the  Beta  Pictoris  disk 
was  likely  in  the  process  of  planetary  formation  due  to  the  presence  of  the  inner  hole,  which 
is  best  explained  through  clearing  mechanisms  associated  with  the  presence  of  planets,  and 
the  fact  that  the  circumstellar  disk  takes  the  form  of  a  flattened  disk,  which  is  consistent 
with  theories  of  planetary  formation. 


Table  5-1:  Index  for  Radial  Brightness  Profile 


Authors 

Range(AU) 

Band 

NE  slope 

SE  slope 

FWHM 

Smith  and  Terrile  (1984) 

120-480 

R 

-4.3 

-4.3 

Parsce  and  Burrows  (1987) 

120-960 

BVRI 

Art3anowicz  et  al.  (1989) 

-3.6 

-3.6 

increasing  with  r 

Lecavalier  des  Etangs  (1993) 

48-360 

BVRI 

-3.6 

-3.6 

Golimowski  et  al.  (1993) 

48-120 

120-360 

R 

-2.38  ±  0.72 
-3.7  ib  0.1 

-1.91  zb  0.89 
-4.0  zb  0.2 

increasing  with  r 

Kalas  and  Jewitt  (1995) 

60-120 

R 

-2.4  rb  0.24 

-2.47  ±  0.36 

43  AU 

Moulliet  et  al.  (1997) 

30-48 

48-77 

77-120 

K’ 

-1.0  ±  0.8 
-1.3  zb  0.2 
-2.5  zb  0.4 

-1.0  zb  0.8 
-1.2  ±  0.2 
-2.9  ±  0.5 

24  AU 

24  AU 

24  AU 

Heap  et  al.  (2000) 

32-74 

56-120 

-1.28  dz  0.04 
-1.79  zb  0.01 

-1.12  zb  0.01 

-1.74  zb  0.04 

^  20  AU 
-  20  AU 

The  above  distances  were  corrected  to  reflect  the  Hipparcos  distance  to  Beta  Pictoris  of  19.28  pc 


The  surface  brightness  profile  measured  by  Smith  and  Terrile  decreased  as  a  power  law 
with  radius  with  an  index  of -4.3  over  the  range  of  120  to  480  AU.  Additional  measurements 
(Parsce  and  Burrows,  1987;  Artymowicz  et  al.,  1989;  Lecavalier  des  Etangs,  1993)  were  made 
that  further  refined  the  power  law  index  for  the  brightness  profile  beyond  120  AU.  The  later 
measurements  showed  an  asymmetry  between  the  extensions  in  both  total  brightness  and 
in  the  power  law  describing  the  brightness  profile.  A  summary  of  the  indices  for  the  power 
laws  describing  the  radial  brightness  distribution  as  determined  in  diflferent  observations  is 
given  in  table  5-1 


58 


Figure  5-3:  Kalas  and  Jewitt  (1995):  (left)  Coronagraphic  image  of  Beta  Pictoris  in  scat¬ 
tered  light,  (right)  Contour  plot  of  surface  brightness  of  the  NE  and  SW  extensions.  The 
’butterfly’  asymmetry  can  be  seen  as  the  asymmetry  in  width  perpendicular  to  the  mid¬ 
plane. 


The  first  images  of  the  disk  within  6  arcsec  (120  AU)  were  produced  using  observations 
made  with  the  Adaptive  Optics  Coronagraph  (AOC)  by  Golimowski  et  al.  (1993).  The  R- 
band  images  showed  a  transition  in  the  power  law  describing  the  surface  brightness  of  the 
mid-plane  at  6  arcsec  (120  AU).  The  change  in  surface  brightness  indicated  a  change  in 
the  natmre  of  the  disk  at  this  distance.  This  change  was  either  in  the  form  of  the  optical 
properties  of  the  material  and/or  a  change  in  the  number  density.  This  conclusion  was 
consistent  with  models  developed  earlier  (Backman  et  al.,  1993)  describing  the  Beta  Pic 
disk  as  having  two  components;  the  inner  disk  (the  region  within  ~  100  AU)  having  lower 
density  and  either  smaller  particles  or  a  less  steep  spatial  gradient. 

The  break  in  the  power  law  describing  the  radial  surface  brightness  distribution  was 
further  verified  by  Kalas  and  Jewitt  (1995)  using  an  i?-band  coronagraph  on  the  University 
of  Hawaii  2.2  m  telescope  on  Mauna  Kea.  Their  image  is  shown  in  figure  5-3  along  with 
a  contour  plot  of  the  surface  brightness.  Once  again,  a  transition  in  the  disk  at  6  arcsec 
was  observed.  Furthermore,  five  major  asymmetries  were  highlighted.  The  first  was  the 
overall  size  of  the  NE  and  SW  extensions  of  the  disk.  The  NE  extension  was  measured  out 
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to  48  arcsec  (925  AU)  and  the  SW  to  only  40  arcsec  (770  AU)  indicating  an  asymmetry  of 
approximately  20%.  The  second  was  an  asymmetry  in  the  overall  brightness  as  measured 
along  the  mid-plane.  At  distances  greater  than  20  arcsec  from  the  star,  the  NE  extension 
is  brighter  along  the  mid-plane  than  is  the  SW  for  the  same  distance  from  the  star.  The 
third  asymmetry  is  in  the  width  of  the  extensions.  At  distances  greater  than  7  arcsec, 
the  FWHM  increases  with  radius.  The  fourth  asymmetry  was  described  as  the  ’butterfly 
asymmetry’  and  referred  to  the  extensions  being  fatter  on  one  side  of  the  mid  plane  than  on 
the  other.  The  fifth  asymmetry  was  described  as  the  ’wing-tilt  asymmetry’  and  referred  to 
the  fact  that  the  mid-planes  of  the  NE  and  SW  extensions  were  inclined  to  one  another  by 
~  1.3  degrees.  Each  of  the  asymmetries  were  scrutinized  to  determine  whether  they  were 
associated  with  a  true  physical  asymmetry  in  the  disk,  or  if  they  represented  an  optical 
effect.  Kalas  and  Jewitt  concluded  that  all  were  true  physical  asymmetries  except  for  the 
wing-tilt  asymmetry  which  was  found  to  be  consistent  with  observing  light  scattered  from  a 
disk  inclined  to  the  line  of  sight.  It  was  also  determined,  based  on  integrated  brightness,  that 
the  overall  number  of  scattering  particles  in  each  extension  were  comparable.  This  implied 
that  the  asymmetries  represented  a  difference  in  the  spatial  distribution  of  dust  rather  than 
a  depletion.  Kalas  and  Jewitt  also  stated  that  the  cause  of  the  asymmetries  remained 
a  mystery.  This  conclusion  was  based  on  considering  two  possible  sources:  a  planetary 
perturbation  or  a  close  stellar  encounter.  The  planetary  scenario  did  not  seem  reasonable 
because  a  planetary  perturbation  was  unlikely  to  cause  asymmetries  on  such  large  scales 
(out  to  ~  1000  AU).  The  stellar  encounter  seemed  unlikely  based  on  a  statistical  analysis 
of  the  probability  associated  with  such  an  event  and  the  fact  that  no  candidate  star  has 
been  observed. 

High  resolution  observations  were  reported  by  Heap  et  al.  (2000)  using  data  obtained 
by  the  Space  Telescope  Imaging  Spectrograh  (STIS)  on  HST.  The  observations  were  made 
using  a  filter  that  covered  200  to  10000  A,  and  achieved  a  resolution  of  ~  0.1”.  The  circum- 
stellar  disk  was  detected  as  close  to  0.75”  which  corresponds  to  15  AU.  The  STIS  images 
showed  a  warped  disk  with  a  maximum  amplitude  of  ^  1.5  AU  at  ^  3.6  arcsec  (70  AU) 
as  seen  in  figure  5-4.  It  was  shown  that  the  observations  could  be  matched  by  modeling 
the  overall  disk  image  as  the  projection  of  two  disks,  one  inclined  to  the  other.  This  was 
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Figure  5-4:  HST  visible  light  image  from  the  Space  Telescope  Imaging  Spectrograph  (STIS) 
clearly  showing  the  warping  of  the  inner  disk  (Heap  et  al.,  2000).  The  bottom  image  shows 
the  same  data  as  the  top,  but  each  vertical  slice  of  the  image  has  been  normalized  to  the 
peak  value  of  that  slice. 
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illustrated  by  decomposing  the  vertical  brightness  profile  of  the  overall  disk  images  into 
two  symmetric  profiles  corresponding  to  a  brighter  main  component  and  a  fainter  inclined 
component.  This  is  illustrated  in  figure  5-5.  The  best  fit  corresponded  to  the  inner,  fainter, 
disk  being  inclined  to  the  main  component  by  4.6  degrees.  Heap  et  al.  also  suggested 
that  the  warp  is  associated  with  the  presence  of  a  planet  and  the  radial  extent  of  the  warp 
depends  on  the  orbital  parameters  of  the  planet  (Larwood  and  Papaloizou,  1997).  The 
analysis  showed  that  the  warp  scaled  with  the  mass  of  the  planet,  the  semi-major  axis  or 
the  planetary  orbit,  and  the  age  of  the  system.  Using  the  estimated  age  of  20  Myr,  the 
planet  ranged  from  17.4  Mj  at  5  AU  to  0.17  Mj  at  50  AU.  The  analysis  was  constrained  by 
the  fact  that  stellar  radial  velocity  variations  are  not  observed  which  implies  Mp/M*  <  0.01. 

5.3.2  Thermal  Emissions 

The  disk  of  Beta  Pictoris  was  first  resolved  in  the  mid-infrared  by  Telesco  et  al.  (1988) 
using  a  20-pixel  bolometer  array  on  the  Infrared  Telescope  Facility  (IRTF)  on  Mauna  Kea. 
The  measurements  made  at  10  and  20  /xm  suggested  sub-micron  particle  sizes  and  an 
inner  hole  with  radius  of  ~  40  AU.  A  later  observation  (Telesco  and  Knacke,  1991)  clearly 
identified  the  presence  of  the  silicate  feature.  Furthermore,  it  was  suggested  that  the  disk 
could  also  include  particles  with  radii  greater  than  1  //m.  An  extension  of  this  work  (Knacke 
et  al.,  1993)  used  extensive  spectroscopic  data  to  constrmn  a  model  of  the  disk.  The  results 
suggested  a  two  component  disk  with  inner  radius  of  30  AU  and  a  transition  region  at  80 
AU.  It  was  also  concluded  that  the  silicates  in  the  disk  consisted  of  a  mixture  of  materials 
similar  to  that  observed  in  comets. 

Lagage  and  Pantin  (1994)  presented  an  infrared  image  of  the  inner  part  of  the  Beta 
Pictoris  dust  disk  by  observing  the  thermal  emissions  at  11.9  fua.  using  the  3.6  m  ESO 
telescope.  The  image  is  shown  in  figure  5-6.  They  achieved  a  spatial  resolution  of  0.3 
arcsec  which  corresponds  to  ~  5  AU.  The  observations  showed  the  dust  was  asymmetrically 
distributed  around  the  stcir  and  depleted  within  2.5  arcsec  (48  AU).  Their  analysis  showed 
the  dust  density  distribution  experienced  a  transition  at  6  arcsec  (120  AU)  which  was 
attributed  to  the  sublimation  of  ice  particles.  They  suggested  that  the  inner  region,  <  48 
AU,  experienced  a  significant  depletion  in  dust,  but  no  strict  void.  This  inner  depletion  was 
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Figure  5-5:  (top)  Decomposition  of  the  vertical  brightness  profile  in  to  two  symmetric 
profiles  associated  with  the  faint  inner  disk  and  the  main  disk,  (bottom)  Plot  illustrating 
the  inclination  of  the  inner  disk  relative  to  the  outer  main  disk  (Heap  et  ah,  2000). 


Figure  5-6:  (left)  11.9  )um  image  of  thermal  emissions  from  Beta  Pictoris  (Lagage  and 
Pantin,  1994).  (right)  2.120  fj,m  coronagraphic  image  using  the  ESO  3.6  m  telescope  with 
adaptive  optics  (Moulliet  et  al.,  1997). 

attributed  to  the  dust  being  trapped  in  gravitational  resonances  or  being  swept  up  by  a 
planet.  They  also  provided  a  qualitative  argument  that  a  planet  with  a  modest  eccentricity 
could  produce  the  asymmetries  observed. 

The  inner  disk  was  later  detected  down  to  1.5  arcsec  from  the  star  and  extending 
out  to  6  arcsec  using  a  coronagraph  at  2.120  /xm  on  the  ESO  3.6  m  with  adaptive  optics 
(Moulliet  et  al.,  1997).  Previous  observations  showed  that  the  radial  surface  brightness 
began  a  transition  at  ~  120  AU,  but  Mouillet  et  al.  showed  that  the  distribution  not 
only  experienced  a  transition  at  120  AU  but  another  at  ~  65  AU  such  that  a  single  power 
law  could  not  reasonable  fit  the  entire  radial  distribution.  The  power  law  developed  to  fit 
the  observations  indicated  that  at  6  arcsec  (120  AU)  the  distribution  becomes  flatter  as 
it  approaches  the  star.  It  was  also  shown  that  the  FWHM  of  the  vertical  profile  stayed 
constant  over  2  arcsec  to  5  arcsec. 

Pantin  et  al.  (1997)  reported  mid-infrared  observations  of  the  inner  disk  which  they 
used  to  derive  the  0-120  AU  density  profile.  The  observations  were  made  using  the  Thermal 
Infrared  Multi  Mode  Instrument  (TIMMI)  mounted  on  the  ESO  3.6  m  telescope  at  11.9 
fim.  The  image  is  shown  in  figure  5-7.  The  density  profile  was  calculated  by  assuming  the 
composition  of  the  dust  material  and  in  turn  developing  a  thermal  model  of  the  disk.  The 


Figure  5-7;  12  /im  image  of  Beta  Pictoris  using  the  Thermal  Infrared  Multi  Mode  Instru¬ 
ment  (TIMMI)  on  ESO  3.6  m  (Pantin  et  al.,  1997). 

material  composition  was  based  on  prior  observations  (Telesco  et  al.,  1988;  Knacke  et  al., 
1993;  Aitken  et  al.,  1993)  mostly  related  to  detection  of  the  silicate  feature.  A  model  based 
on  a  single  grain  composition  proved  unable  to  make  a  good  fit  to  the  data,  but  a  model  with 
two  components  was  more  successful.  The  transition  between  the  components  was  defined 
by  the  ice  boundary  estimated  to  be  at  ~  90  AU  beyond  which  the  particles  contain  some  ice. 

5.3.3  Spectroscopic  Observations 

An  effort  to  combine  ground  and  space  based  spectroscopic  observations  of  the  Beta 
Pictoris  disk  was  performed  by  Lagrange  et  al.  (1996)  to  test  the  Falling  Evaporating  Body 
(FEB)  scenario.  The  FEB  theory  proposes  that  comet  sized  bodies  on  near  circular  orbits 
become  star-grazers  due  to  planetary  perturbations.  These  bodies  are  generally  thought  to 
be  asteroid-like  with  an  icy  nucleus  rather  than  fully  icy  comets  (Karmann  et  al.,  2001). 
As  these  bodies  move  near  the  star  they  begin  to  evaporate.  If  this  occms  during  a  transit 
across  the  line  of  sight,  it  will  produce  observable  absorption  lines  with  short  timescale 
variations. 

Velocities  and  timescales  for  variability  were  determined  by  looking  at  absorption  lines 
associated  with  Ca  II,  Mn  II,  Fe  II,  Al  III,  and  Al  II.  Their  approach  was  based  on  correlating 
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the  observed  red-shift  with  variability  timescales,  line  width,  and  line  depth.  Two  regimes 
were  reported:  1)  low  velocity,  slowly  varying  lines  and  2)  high  velocity,  rapidly  varying 
lines.  The  high  velocity  lines  were  determined  to  be  consistent  with  and  well  explained  by 
the  evaporation  of  single  bodies  within  10  stellar  radii.  The  low  velocity  lines  were  best 
explained  by  the  evaporation  of  families  of  bodies  crossing  the  line  of  sight  at  large  distances, 
~  20  -  30  stellar  radii,  from  the  star.  Therefore  the  conclusion  reached  by  Lagrange  et  al. 
was  that  the  variable  features  could  be  interpreted  as  signs  of  km-sized  solid  bodies  on  star 
grazing  orbits  evaporating  near  periastron  -  the  FEB  scenario.  It  was  also  argued  that, 
because  all  of  the  variability  was  seen  on  relatively  short  time  scales,  mean  motion  and  not 
secular  resonances  were  the  dominant  factor  (Lagrange  et  al.,  1996). 

It  was  also  noted  that  the  ratio  between  red-shifted  and  blue-shifted  events  was  not  1:1. 
In  fact,  for  every  one  blue-shifted  event  observed,  nine  red-shifted  events  were  recorded.  It 
was  suggested  (Artymowicz,  1997)  that  this  was  indicative  of  a  preferential  orientation  for 
the  orbits  of  the  solid  bodies. 

5.3.4  Photometric  Variations 

Photometric  variations  of  Beta  Pictoris  were  studied  (Lecavaher  des  Etangs  et  al.,  1995) 
using  data  compiled  from  near-continuous  observations  made  between  1975  and  1992.  A 
variation  was  detected  with  a  timescales  of  ~  1000  days  and  an  amplitudes  of  ~  0.01 
magnitudes.  It  was  suggested  that  this  variation  was  attributed  to  inhomogeneities  in  the 
azimuthal  distribution  of  dust  possibly  caused  by  the  presence  of  a  planet  or  planetesimals. 
Another  short  term  variation  was  observed  in  November  of  1981.  The  timescale  associated 
with  this  variation  was  ~  30  days  and  reached  a  maximum  amplitude  of  0.04  magnitudes. 
It  was  proposed  that  this  short  time  scale  variation  was  the  result  of  a  passage  of  a  planet  - 
an  occultation  of  Beta  Pictoris  by  a  planet  roughly  the  size  of  Jupiter.  The  region  around 
the  planet  was  considered  to  be  relatively  clear  of  dust  due  to  accretion  or  ejection  of  the 
particles  and  this  resulted  in  the  brightness  enhancement. 
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5.4  Dust  Properties 

This  section  provides  a  brief  review  of  the  properties  of  the  dust  comprising  the  Beta 
Pictoris  circumstellar  disk.  It  is  these  properties,  based  on  observation,  that  will  be  used 
in  the  models  described  in  the  following  chapters. 

5.4.1  Dust  optical  properties 

The  presence  of  the  silicate  feature  first  identified  by  Telesco  and  Knacke  (1991)  in¬ 
dicated  that  the  particles  were  likely  similar  to  material  common  within  the  solar  system. 
Prior  to  1998  the  material  was  generally  modeled  as  Mie  spheres  of  astronomical  silicates 
whose  optical  material  properties  were  described  by  Draine  and  Lee  (1984).  However,  Li 
and  Greenberg  (1998)  performed  a  more  detailed  analysis  of  the  optical  properties  of  the 
dust  within  the  Beta  Pictoris  system  and  concluded  that  the  dust  had  a  more  complicated 
structure  and  hence  interacted  with  light  in  a  way  much  different  than  what  would  be 
expected  for  homogeneous  spheres. 

Li  and  Greenberg  (1998)  proposed  that  the  dust  material  within  Beta  Pictoris  is 
cometary  in  origin  and  consisted  of  aggregates  of  particles  with  silicate  cores  and  an  organic 
refractory  mantle.  Furthermore,  these  particles  are  highly  porous,  ~  P  >  0.95.  It  was  also 
suggested  that  in  the  outer  regions,  the  particles  possessed  a  water  dominated  ice-mantle. 
In  general  Li  and  Greenberg  asserted  that  the  material  is  created  near  the  central  star  due 
to  the  evaporation  of  the  comets  and  distributed  throughout  the  disk  via  radiation  pressure. 

5.4.2  Dust  Particle  Size 

Typical  particle  sizes  determined  firom  past  observations  and  analysis  are  thought  to 
have  a  radius  of  «  1  to  3  /jim  (Artymowicz,  1997).  This  estimate  weus  beised  on  both  dynam¬ 
ical  considerations  and  constraints  firom  the  10  fim  emissions.  Artymowicz  suggested  that 
particles  smaller  than  ~  1  /xm  are  quickly  removed  from  the  system  by  radiation  pressure. 
Furthermore,  the  particles  must  be  ~2  /xm  to  stay  in  orbit  long  enough  to  produce  the  10 
//m  emissions.  Further  analysis  of  the  typical  particle  size  is  addressed  in  chapter  7. 


Figure  5-8:  HST  WFPC2  and  UH  images  after  subtraction  of  a  smooth  model  disk.  The 
lines  show  the  corresponding  brightness  enhancements  located  at  ~  750  AU  (Kalas  et  ah, 
2000). 


5.5  Models 

The  most  extensive  numerical  modeling  done  on  the  dust  disk  of  Beta  Pictoris  has 
been  based  on  the  presence  of  a  single  planet  inclined  to  the  plane  of  the  disk  or  due  to 
the  close  encounter  with  another  star.  An  alternative  approach  was  based  on  considering 
a  driving  force  produced  by  stellar  radiation  (Armitage  and  Pringle,  1996).  In  this  theory, 
stars  with  luminosities  greater  than  lOL©  are  expected  to  produce  warps  in  circumstellar 
disks.  However,  when  the  Hipparcos  data  determined  a  new  distance  to  Beta  Pictoris,  and 
thus  a  new  value  for  luminosity,  this  theory  no  longer  seemed  applicable  (Heap  et  al.,  2000). 

This  section  discusses  first  the  theory  of  a  close  stellar  encounter  and  then  the  models 
based  on  a  single  planet  on  an  inclined  orbit. 

5.5.1  Close  Stellar  Encounter 

Kalas  et  al.  (2000)  compiled  images  from  three  different  telescopes  using  four  different 
instruments:  Las  Campanas  2.5  m.  University  of  Hawaii  2.2  m,  HST  Planetary  Camera, 
and  HST  Wide  Field  Camera.  After  subtracting  a  smooth,  symmetric  model  disk,  bright 
features  located  ~  3.9  arcsec  (~  750)  AU  from  the  star  in  the  NE  extension  were  identified  in 
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Figure  5-9:  Results  of  a  numerical  simulation  of  a  close  stellar  encounter  with  a  disk  of 
planetesimals.  (Kalas  et  al.,  2000). 

all  four  images.  These  brightness  enhancements  were  interpreted  as  signs  of  a  ring  system. 
Furthermore,  since  the  features  were  not  seen  in  the  SW  extension,  it  was  proposed  that 
the  rings  were  not  centered  on  the  star.  Therefore,  it  was  suggested  that  these  features 
were  the  edge-on  view  of  an  asymmetric  ring  system  within  the  disk.  Two  of  the  images 
used  for  the  analysis  are  shown  in  figure  5-8. 

To  explain  the  nature  of  the  ring  system,  numerical  simulations  were  performed  to 
study  the  effect  of  a  close  stellar  encounter  with  a  quiescent  disk  of  planetesimals  (Kalas  et 
al.,  2000;  Larwood  and  Kalas,  2001).  The  model  included  10®  collisionless  particles  initially 
in  circular  orbits  about  a  point-mass  potential.  The  stellar  flyby  consisted  of  a  secondary 
point  mass,  M  ~  0.3M*,  on  a  parabolic  orbit  which  was  not  coplanar  to  the  original  plane 
of  the  disk. 

The  encounter  resulted  in  the  planetesimals  being  scattered  outward  several  times  the 
initial  disk  radius  but  maintaining  a  well  defined  disk  structmre.  Large  scale  ring  structures 
formed  but  dissolved  after  about  10®  years,  or  10  orbital  periods  at  500  AU,  due  to  phase 
mixing  -  the  dispersal  of  particles  based  on  differences  in  their  periodicity.  It  was  suggested 
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that  the  observed  ring  structure  was  a  tightly  wound  spiral  pattern  resulting  from  the 
induced  motion  of  the  primary  star  by  the  perturber  (Larwood  and  Kalas,  2001).  The 
results  of  one  simulation  are  shown  in  figure  5-9. 

Multiple  simulations  were  performed  and  artificial  isophotes  were  created  and  compared 
to  observations  of  Beta  Pic.  The  best  fit  occurred  for  an  encounter  distance  of  ~  700  AU. 
For  this  value  the  probability  of  such  an  encounter  is  very  low,  ~  0.01%  in  10®  years  (Kalas 
et  ah,  2000).  Furthermore,  no  suitable  candidate  has  been  identified.  Kalas  et  al.  (2000) 
suggested  that  if  Beta  Pic  formed  with  a  bound  companion  that  the  probability  increased 
by  am  order  of  magnitude.  However,  no  companion  for  Beta  Pictoris  has  been  identified. 

Kalas  et  ad.  (2000)  ailso  suggested  that  a  close  stellar  encounter  could  provide  an  ex¬ 
planation  for  the  observed  warp  in  the  disk.  It  was  argued  that  a  non-coplanar  flyby  would 
increase  the  inclination  and  eccentricity  of  the  planetesimal  orbits.  This  would  create  a 
second  plane  whose  edge  on  projection,  when  coupled  with  the  main  disk,  would  create  the 
appearance  of  a  warp. 

5.5.2  Single  Planet  on  an  Inclined  Orbit 

Based  on  the  observation  of  a  possible  occultation  of  Beta  Pictoris  (Lecavalier  des 
Etangs  et  al.,  1995)  and  the  observed  high  rate  of  cometary  infall  (Lagrange  et  al.,  1996), 
Moulliet  et  al.  (1997b)  suggested  that  a  planet  was  likely  present  and  performed  numerical 
simulations  to  illustrate  the  effects  of  a  single  planet  on  a  circumstellar  disk.  The  model 
was  based  on  smoothed  particle  hydrodynamics  with  pressure  and  viscous  forces  removed. 
They  also  made  the  assumption  that  the  disc  shape,  observed  through  the  small  particles, 
was  driven  by  gravitational  effects  on  larger,  longer  lived  particles. 

The  disk  model  consisted  of  15,000  particles  initially  on  circular  orbits  populated  in 
such  a  manner  as  to  be  consistent  with  an  assumed  matter  distribution.  A  planet  was 
introduced  with  an  eccentric  orbit  inclined  to  the  disk  mid-plane.  Moulliet  et  al.  (1997) 
suggested  that  as  a  result  of  the  break  in  symmetry,  a  deformation  in  the  disk  occurs  that 
is  directly  related  to  the  planets  inclination.  The  results  of  one  simulation  are  shown  in 
figure  5-10. 
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Figure  5-10:  Results  of  a  numerical  model  involving  15,000  particles  in  a  disk  perturbed 
by  a  planet  on  an  inclined  orbit  (Moulliet  et  al.,  1997). 


By  calculating  the  flux  due  to  stellar  light  being  scattered  from  the  dust  in  the  disk 
and  integrating  along  a  line  of  sight,  Moulliet  et  al.  (1997b)  determined  an  apparent  sur¬ 
face  brightness  distribution  which  they  compared  to  previous  coronagraphic  observations 
(Moulliet  et  al.,  1997).  Using  the  constraints  associated  with  the  observations  and  the  fact 
that  no  radial  velocity  variations  have  been  observed  for  Beta  Pictoris,  it  was  determined 
that  a  planet  should  be  located  between  1  and  20  AU,  be  inclined  at  3  -  5°,  and  have  a  mass 
between  10"^  and  10“®M*.  It  was  also  suggested  that  the  planet  should  have  a  significant 
eccentricity  to  explain  the  high  rate  of  cometary  infalls. 

Augereau  et  al.  (2001)  presented  a  ’complete  model’  of  the  Beta  Pictoris  disk  that 
successfully  accounted  for  the  major  asymmetries  seen  in  both  scattered  light  and  thermal 
emissions.  The  model  was  based  on  a  disk  of  planetesimals  extending  out  to  150  AU  per¬ 
turbed  by  a  single  planet  on  an  inclined  orbit.  The  small  scale  asymmetries,  <  120  AU, 
were  explained  through  the  gravitational  perturbation  of  the  planet.  The  large  scale  asym¬ 
metries  were  explained  by  radiation  pressure  extending  the  distribution  of  small  particles, 
created  by  collisional  processes,  far  beyond  that  of  the  parent  body  disk. 

The  modeling  approach  was  to  construct  a  parent  disk  from  100  concentric  rings  each 
composed  of  10^  particles.  The  disk  was  allowed  to  phase  mix  for  20  orbital  periods.  The 


Figure  5-11:  A  plot  of  the  vertical  position  of  the  mean  parent  body  disk  mass  (Augereau 
et  al.,  2001).  The  maximum  is  at  a  distance  corresponding  to  ~  75  AU. 

orbital  planes  were  precessed  based  on  the  local  precession  rate  determined  by  the  dynamical 
influence  of  the  planet.  The  disk  W2is  then  allowed  to  phase  mix  again.  It  was  then  assumed 
that  each  parent  body  produced  a  distribution  of  grains  each  with  a  speciflc  value  of  the 
ratio  between  the  force  due  to  radiation  pressure  and  the  force  due  to  gravity  from  the 
star.  The  spatial  distribution  of  small  particles  was  determined  for  11  difierent  values  of 
running  between  0.001  and  0.45.  Simulated  images  were  made  for  each  value  of  /5.  The 
final  image  was  a  superposition  using  weighting  factors  based  on  the  assumed  grain  size 
distribution. 

Augereau  et  al.  (2001)  claimed  that  the  presence  of  the  planet  causes  the  parent  body 
disk  to  process  differentially.  The  inner  region  of  the  disk  becomes  coplanar  with  the  planet 
due  to  its  more  rapid  rate  of  precession.  This  means  that  the  extent  of  the  warp  is  based 
on  the  processional  rate  and  the  age  of  the  system,  suggesting  that  the  warp  propagates 
outward.  A  plot  illustrating  how  the  mass  of  the  parent  body  disk  was  distributed  is  shown 
figure  5-11. 

The  best  fit  to  observations  occurred  with  a  perturbing  planet  of  mass  M/M*  =  10“^, 
semi  major  axis,  o  =  10  AU,  e  =  0.0,  and  7  =  3°. 

5.5.3  Limitations  of  Models 

The  modeling  of  the  Beta  Pictoris  disk  that  has  been  discussed  in  this  chapter  has  been 
based  on  a  single  planet  with  an  inclination  relative  to  the  disk.  However,  these  models 
do  not  address  the  issue  of  how  such  a  configuration  could  airise.  An  alternative  approach 
is  to  consider  a  system  of  two  or  more  planets  that  could  produce  secular  perturbations. 


Figure  5-12:  Visualization  of  a  SIMUL  model  representing  the  distribution  of  surface  area. 
The  upper  left  frame  shows  the  bottom  of  the  disk,  with  the  disk  tilted  30  degrees  from 
edge-on.  Going  left  to  right  and  top  to  bottom,  the  tilt  angle  -15,  -3,  0,  15,  20,  30,  45,  and 
60  degrees.  The  images  clearly  show  the  resulting  two  disk  structure. 

As  was  discussed  in  chapter  3,  secular  perturbation  theory  is  a  well-developed  theory  used 
for  analyzing  the  long-term  orbital  perturbations  in  a  many-body  system.  This  theoretical 
framework  has  been  applied  extensively  to  problems  within  the  solar  system.  In  fact, 
analysis  of  the  secular  perturbations  arising  from  Jupiter  and  Saturn  can  be  shown  to 
create  a  forced  inclination  that  varies  with  semi-major  axis,  o,  (see  for  example  Dermott  et 
al.  (1984,  1985);  Mmray  and  Dermott  (1999)).  It  can  also  be  shown  that  this  variation  of 
forced  inclination  can  produce  a  warp  in  a  dust  disk. 

«r 

In  this  chapter  it  was  shown  that  warps  can  be  simulated  by  considering  the  edge-on 
view  of  a  multi  component  system  in  which  the  components  are  concentric  and  have  a  non¬ 
zero  mutual  inclination.  Therefore,  any  modeling  technique  utilized  must  produce  results 
consistent  with  this  effect.  In  other  words,  the  approach  used  in  this  research  should  ideally 
reproduce  the  multi-component  system. 
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One  of  the  tools  utilized  in  this  research,  which  will  be  discussed  in  detail  later,  calcu¬ 
lates  the  distribution  of  surface  area  based  on  the  functional  form  of  the  orbital  parameters 
with  respect  to  semi-major  axis.  One  such  realization  is  shown  in  figure  5-12.  This  image 
was  created  using  Beta  Pictoris  with  a  two  planet  system.  The  orbital  parameters  were 
those  of  Jupiter  and  Saturn,  but  with  the  Saturn-like  planet  pushed  out  to  20  AU.  This 
image  illustrates  how  the  planes  of  symmetry  of  a  disk  model  are  infiuenced  by  secular  per¬ 
turbations  associated  with  planets  of  relatively  low  inclination.  Furthermore,  it  illustrates 
how  such  a  configuration  can  produce  a  structure  consisting  of  two  disks,  one  inclined  to  the 
other.  These  results  are  the  same  basic  results  achieved  in  the  other  models,  but  without 
the  less  plausible  scenario  of  having  a  single  planet  on  an  orbit  inchned  the  disk  from  which 
it  formed. 


CHAPTER  6 

THE  MODELING  PROCESS 


6.1  Introduction 

Previous  modeling  of  the  Beta  Pictoris  dust  disk  consisted  primarily  of  N-body  or  hy¬ 
drodynamic  analysis.  These  approaches  consist  of  assuming  an  initial  configuration  for  the 
dust  particles  or  their  parent  bodies,  simulating  key  physical  processes,  aad  allowing  the 
system  to  evolve  with  time  from  to  to  t.  The  approach  taken  in  this  research  is  different  in 
that  the  simulation  process  does  not  track  and  record  the  evolution  of  the  particle  orbits 
or  the  distribution  of  particles  within  a  region  of  space,  but  instead  attempts  to  directly 
calculate  the  overall  effect  of  gravitational  perturbations  on  long  timescales.  In  other  words, 
the  modeling  described  here  directly  calculates  the  configuration  at  time  t.  This  chapter 
will  discuss  the  simulation  process  utilized  in  this  research  in  developing  a  disk  model. 

6.2  SIMUL 

Dynamic  modeling  of  the  dust  cloud  within  the  solar  system,  the  zodiacal  cloud,  has 
been  studied  for  many  years  and  many  tools  have  been  developed  to  facilitate  this  research. 
In  particular,  a  program  called  SIMUL  was  developed  initially  by  Dermott  and  Nicholson 
(1989)  and  then  by  Xu  et  al.  (1993)  to  construct  a  three-dimensional  numerical  model  of  the 
solar  system  dust  bands  and  illustrate  their  link  to  the  major  Hirayama  asteroid  families. 

Key  to  SIMUL  is  the  idea  that  a  disk  is  a  distribution  of  many  dust  particle  orbits,  along 
which  the  dust  particles  are  distributed.  The  characteristics  of  an  orbit  are  defined  and  or 
influenced  by  both  initial  conditions  and  other  perturbing  forces.  The  dust  particle  orbits 
are  described  by  the  standard  orbital  elements  of  semi  major  axis,  eccentricity,  longitude 
of  pericenter,  inclination,  and  longitude  of  ascending  node.  Furthermore,  the  observable 
orbital  elements,  the  osculating  elements,  are  associated  with  two  components:  the  forced 
elements  and  the  proper  elements. 
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As  stated  before,  the  forced  and  proper  elements  can  be  linked  to  the  perturbing 
planets  or  the  conditions  in  which  the  dust  was  created.  This  implies  that  characteristics  of 
a  particle’s  orbit  can  be  determined  by  knowing  how  the  dust  composing  the  disk  was  created 
and  how  the  orbit  is  perturbed  by  massive  bodies.  Furthermore,  these  two  components  are 
not  coupled  and  can  be  determined  independently.  The  forced  elements  are  determined  by 
performing  the  calculations  described  in  chapter  3.  The  proper  elements  are  determined 
by  making  assumptions  based  on  what  we  know  of  dust  creation  within  the  solar  system 
extrapolated  to  other  environments. 

The  distribution  of  the  orbits  in  space  is  determined  by  both  the  distribution  of  the 
bodies  creating  the  dust  and  the  way  in  which  the  dust  is  affected  by  radiation  affects. 
The  dust  is  generally  assumed  to  be  created  from  the  destruction  of  larger  bodies.  After 
creation,  the  dust,  depending  on  the  grain  size  and  density,  will  either  stay  in  the  same 
orbit,  be  blown  out  of  the  system  by  radiation  pressure,  or  spiral  into  the  star  by  PR  drag. 
These  factors  will  determine  the  spatial  distribution  of  dust  of  a  given  size  as  a  function  of 
distance  from  the  star. 

The  primary  function  of  SIMUL  is  to  produce  a  three  dimensional  mapping  of  total 
surface  area  associated  with  dust  particles  for  any  specified  value  of  Beta.  As  input,  SIMUL 
requires  either  tabulated  or  functional  forms  of  forced  and  proper  orbital  elements.  It  also 
requires  the  radial  distribution  of  dust  particle  orbits  which  is  influenced  by  the  interaction 
of  the  dust  with  the  radiation  from  the  star. 

The  models  produced  by  SIMUL  are  based  on  dividing  up  the  total  surface  area,  a 
model  parameter,  among  a  large  number  of  dust  particle  orbits.  For  each  orbit,  the  surface 
area  is  distributed  evenly  in  mean  anomaly.  Orbits  are  determined  by  randomly  selecting 
a  value  for  semi-major  axis  based  on  a  distribution  function  determined  from  the  expected 
distribution  of  dust  particle  orbits.  For  the  specific  value  of  semi-major  axis  selected,  the 
corresponding  forced  and  proper  elements  are  selected  from  input  data,  and  the  orbit  is 
mapped  out  in  three  dimensional  space.  The  orbit  thus  defined,  contains  all  characteristics 
associated  with  the  initial  creation  of  the  dust  particle  and  the  effects  of  perturbations  over 
long  timescales.  This  process  is  repeated  any  number  of  times,  ~  10®  in  this  research,  until 
the  disk  is  adequately  populated  with  enough  orbits  to  reproduce  the  specified  total  surface 


76 


area.  The  final  spatial  distribution  is  generated  by  determining  where  each  of  the  orbits 
intersects  an  element  of  a  predefined  three  dimensional  grid. 

6.3  Visualization 

The  final  output  of  SIMUL  is  a  three  dimensional  array  in  which  each  element  repre¬ 
sents  the  surface  area  per  volume  for  a  given  particle  diameter.  To  transform  this  result  into 
a  simulated  dust  distribution  several  steps  must  be  taken.  If  a  mapping  of  the  distribution 
of  surface  area  is  desired,  then  integration  along  a  line  of  sight,  through  the  disk  described 
by  the  three  dimensional  array  output  by  SIMUL  can  be  done  directly  after  the  appropriate 
rotations  are  performed.  If  a  model  disk  is  desired,  then  further  input  is  required. 

An  actual  circumstellar  disk  will  contain  dust  that  is  emitting  thermal  radiation.  This 
dust  is  heated,  primarily  by  the  short  wavelength  radiation  of  the  star,  and  emits  in  the 
long  wavelength  region  of  the  spectrum.  This  makes  the  mid-IR  ideal  for  observing  the 
dust  in  circumstellar  disks.  If  a  simulation  of  these  emissions  is  desired,  the  thermal  and 
optical  properties  of  the  dust  must  be  considered.  The  properties  that  determine  how  a 
particle  absorbs  light,  is  heated,  and  then  re-radiates  energy  is  based  on  the  dust’s  material 
properties  and  structure.  Coupling  assumptions  about  these  properties  with  the  properties 
of  the  star  allow  estimates  to  be  made  about  the  temperature  of  dust  particles  as  a  function 
of  distance  from  the  star  and  particle  size.  This  in  turn  produces  an  estimated  flux  per  unit 
of  surface  area  for  a  particular  sized  particle. 

Finally,  this  flux  can  be  utilized  with  the  surface  area  model  to  determine  an  overall 
flux  distribution  for  the  entire  disk.  This  is  done  by  selecting  an  observing  direction  with  re¬ 
spect  to  the  disk,  and  integrating  along  that  line-of-sight.  This  produces  a  two-dimensional 
mapping  for  which  each  pixel  value  is  related  to  the  total  flux  from  that  location. 

6.4  Summary 

The  approach  to  using  secular  perturbation  theory  to  model  the  warp  in  a  disk  like 
that  of  Beta  Pictoris  requires  determining  the  properties  of  a  planetary  system  that  will 
produce  a  disk  model  that  has  features  matching  those  seen  in  the  observed  disk.  This 
requires  first  choosing  a  configuration  of  planets  and  calculating  the  forced  orbital  elements 
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as  a  function  of  semi-major  axis.  Then,  by  considering  the  other  physical  processes  involved 
in  determining  the  distribution  of  material  in  the  disk,  a  model  disk  can  developed.  The 
following  chapter  illustrates  how  the  characteristics  and  features  quantified  in  an  observation 
can  constrain  the  input  parameters  in  this  process  and  produce  an  output  that  matches  the 
original  observation. 

A  subtle  but  important  point  to  emphasize  about  this  approach  is  that  the  dynamics 
of  the  dust  particles  themselves  are  not  being  modeled.  It  is  actually  the  larger  and  more 
massive  source  bodies  that  axe  being  modeled  dynamically.  This  is  necessary  as  the  parti¬ 
cles  that  are  generally  responsible  for  the  majority  of  the  thermal  emissions  are  measured 
on  the  scale  of  microns.  These  particles  are  generally  short  lived  in  a  system  like  Beta 
Pictoris,  as  they  are  quickly  expelled  firom  the  system  on  parabolic  or  hyperbolic  orbits 
due  to  radiation  pressure.  However,  as  will  be  discussed  in  the  following  chapter,  under 
the  appropriate  conditions  it  is  valid  to  eissume  that  the  observed  particles  trace  the  larger 
unseen  bodies.  This  is  because  the  spatial  distribution  of  the  small  particles,  despite  being 
quickly  blown  out,  share  many  of  the  characteristics  of  the  spatial  distribution  of  the  source 
bodies  from  which  they  were  created.  This  is  discussed  further  in  the  next  chapter. 


CHAPTER  7 

CHARACTERIZING  THE 
BETA  PICTORIS  DISK 

7.1  Introduction 

The  recent  mid-IR  images  of  Beta  Pictoris  show  with  great  detail  the  features  and 
asymmetries  of  a  warped  circumstellar  disk.  Such  images  allow  the  unique  opportunity  to 
directly  observe  the  effects  of  dynamical  perturbations  on  the  disk  material.  This  resemch 
asserts  that  these  unique  features  and  asymmetries  are  linked  to  the  presence  of  a  multiple 
planet  system. 

The  goal  of  the  next  two  chapters  is  to  illustrate  how  the  theory  and  modeling  process 
presented  in  previous  chapters  can  be  applied  in  such  a  way  as  to  reproduce  the  important 
characteristics  of  a  warped  circumstellar  disk  such  as  those  determined  from  the  observations 
of  Beta  Pictoris.  Considering  the  fact  that  only  a  single  observing  perspective  is  available,  a 
unique  solution  cannot  be  found  with  the  number  of  free  parameters  in  a  system  as  complex 
as  that  of  Beta  Pictoris.  Nevertheless,  the  next  two  chapters  will  develop  a  logical  approach 
to  constrain  a  proposed  system.  Furthermore,  it  will  be  shown  that  at  least  one  plausible 
configuration  of  planets  will  produce  the  same  t3rpe  of  interesting  features  that  have  been 
identified  in  several  observations,  thus  adding  validity  to  the  general  approach. 

In  this  chapter,  the  recent  observations  of  Beta  Pictoris  made  in  the  mid-IR  by  the 
University  of  Florida  team  (Telesco  et  al.,  2003)  will  be  discussed.  Similar  observations 
made  by  other  groups  will  also  be  addressed.  Highlighted  in  this  discussion  will  be  the 
identification  of  the  major  symmetry  planes  and  the  associated  relative  inclinations.  The 
major  symmetry  planes  are  the  most  well  studied  feature  of  the  disk  and  will  be  the  first 
characterized  in  the  modeling. 

The  physical  conditions  in  the  circumstellar  environment  place  important  constraints 
on  the  modeling  process.  The  assumptions  made  about  the  physical  processes  specific  to 
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the  environment  around  Beta  Pictoris  are  also  discussed  in  this  chapter.  This  discussion 
includes  analysis  on  the  timescales  of  the  various  physical  processes  that  occur  within  the 
disk.  The  link  between  these  processes  and  environmental  properties  such  as  dust  grain 
composition  and  the  estimated  optical  depth,  is  also  addressed.  Finally,  the  assumptions 
about  the  values  of  the  proper  orbital  elements  will  be  discussed.  In  chapter  4  the  link 
between  the  observable  characteristics  of  the  disk  and  the  proper  orbital  elements  was 
discussed.  This  chapter  will  discuss  the  ways  in  which  the  observation  of  Beta  Pictoris 
constrain  the  values  for  the  proper  elements  that  are  used  in  the  models  to  be  presented  in 
chapter  8.  These  discussions  will  be  used  to  argue  the  validity  of  the  model  assumptions 
and  the  applicability  of  secular  perturbation  theory  to  this  problem. 

7.2  Recent  Mid-IR  Observations  of  Beta  Pictoris 

The  University  of  Florida’s  IR  group  made  10.8  and  18.2  micron  observations  of  the 
Beta  Pictoris  disk  using  the  instrument  OSCIR  on  the  8  m  Gemini  South  telescope  on  Cerro 
Pachon  in  Chile  (Telesco  et  ah,  2003).  The  observations  made  in  December  of  2001  are 
shown  in  figure  7-1.  The  observations  shown  in  this  image  represent  the  highest  resolution 
images  (approximately  0.5  arcsec  at  18  |im)  ever  taken  at  these  mid-IR  wavelengths. 

The  images  expose  emission  from  the  disk  extending  out  to  a  radius  of  approximately 
100  AU  fi-om  the  central  star.  The  images  also  reveal  several  characteristics  that  hint  at 
the  presence  of  a  planetary  system.  Clearly  seen  in  the  18  /xm  image  is  an  ”S-curve”  shape 
in  the  southwest  extension.  This  is  consistent  with  the  warping  seen  in  other  observations. 
The  UF  images  however  indicate  the  details  of  the  warping  in  such  a  way  as  to  expose 
the  complex  shape  of  the  disk’s  mid-plane.  This  hints  at  a  complex  mechanism  being 
responsible  for  shaping  the  spatial  distribution  of  circumstellar  material. 

Also  observed  in  the  images  are  regions  in  the  disk’s  mid-plane  that  suggest  depletions 
or  enhancements  in  the  dust  number  density.  Consistent  with  theories  of  gap  creation  by 
large  bodies  in  a  circumstellar  disk,  these  voids  or  gaps  will  be  used  later  to  constrain  the 
locations  of  the  planets  used  in  the  model.  Identification  of  these  features  is  discussed  later 
in  this  chapter. 


Boto.  Piotoris  lO  miorojas 


0  50  100  150  200 


Asimnomtcal  Uni  b 


Piotoiris  18  i33.ioirow.s 


0  50  100  150  200 

A^iroxtomicol  Unit* 


Figure  7-1;  Mid-IR  images  of  Beta  Pictoris  in  10  fim  (top)  and  18  /im  (bottom). 
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Figure  7-2:  Mid-IR  images  of  Beta  Pictoris  in  12  /xm  (upper)  and  18  p  m  (lower)  (Wein¬ 
berger  et  al.,  2003).  The  images  show  the  inclined  plane  of  the  inner  disk  (^  1  arcsec)  and 
outer  disk  (~  3  arcsec). 
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Weinberger  et  al.  (2003)  and  Wahhaj  et  al.  (2003)  also  presented  mid-IR  observations 
of  Beta  Pictoris  at  11.7  and  17.9  pm  and  17.9  pm  respectively.  Both  observations  were 
made  using  the  Keck  observatory  in  Hawaii.  These  images  illustrate  the  changes  in  relative 
inclination  angle  along  the  mid-plane.  Weinberger  et  al.  identified  the  inner  warp  in  the 
disk  (<  20  AU)  as  well  as  the  previously  identified  large  scale  symmetry  planes.  Wahhaj  et 
al.  not  only  observed  the  warping  of  the  mid-plane,  but  they  also  identified  emission  peaks 
within  the  disk  that  they  associated  with  a  ring  system.  The  estimated  ring  radii  were 
14  ±  1,  28  ±  3,  52  ±  2,  and  82  ±  2  AU.  The  rings  were  also  argued  to  have  inclinations 
that  alternate  in  orientation  relative  to  the  invariant  plane  of  the  disk.  The  Wahhaj  et  al. 
image  is  consistent  with  the  UF  images  in  the  sense  that  it  shows  a  disk  consisting  of  a 
ring-like  structure  but  it  does  show  inconsistencies  in  the  locations  and  inclinations  of  the 
rings.  However,  inconsistencies  should  not  be  entirely  unexpected  as  the  UF  observations 
were  made  at  significantly  lower  airmass  than  that  of  the  Wahhaj  observation. 

Further  analysis  of  the  UF  Beta  Pictoris  images  reveals  essentially  the  same  symmetry 
planes  as  seen  in  the  images  of  Heap  et  al.  (2000),  Weinberger  et  al.  (2003),  and  Wahhaj  et 
al.  (2003).  Together  these  observations  create  a  consistent  picture  of  a  distribution  of  dust 
with  three  major  symmetry  planes.  Furthermore,  the  secondary  effects  of  planets  that  may 
be  responsible  for  the  large  scale  warping  are  seen  within  the  UF  images.  The  modeling  of 
all  of  these  features  is  saved  for  the  next  chapter. 

7.3  The  Physical  Enviromnent 

Modeling  the  Beta  Pictoris  system  in  the  context  of  secular  perturbation  theory  re¬ 
quired  three  primary  steps.  First,  the  timescales  associated  with  the  various  processes  were 
determined  to  ensure  that  the  system  has  had  enough  time  to  fully  evolve;  i.e.  all  regions  of 
the  observed  disk  have  been  influenced  by  the  secular  perturbations.  Key  to  understanding 
the  timescales  associated  with  the  disk,  is  a  reasonable  estimate  of  the  dust  properties  - 
composition  and  typical  size.  Second,  the  location  of  the  planets  had  to  be  assumed  using 
characteristics  of  the  observed  morphology.  Finally,  the  forced  elements  needed  to  be  cal¬ 
culated  and  used  to  simulate  a  distribution  of  dust.  This  section  discusses  the  first  step  of 
determining  the  physical  conditions  existing  within  the  disk,  starting  with  an  analysis  of 
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Figure  7-3:  Mid-IR  images  of  Beta  Pictoris  at  17.9  (Wahhaj  et  al.,  2003).  The  image 
shows  the  inclined  plane  of  the  inner  disk  (^20  AU  and  outer  disk  (~  60  AU). 

the  timescales  associated  with  the  processes  that  most  affect  the  morphology  of  the  disk. 

7.3.1  Typical  Particle  Size 

The  interaction  between  radiation  from  the  star  and  the  dust  particles  play  an  impor¬ 
tant  role  in  determining  the  distribution  of  dust  in  a  disk.  For  a  thorough  discussion  on 
the  dust  particle  size  distribution  that  is  expected  for  a  circumstellar  disk  see  Wyatt  et  al 
(1999)  and  references  within.  As  a  general  summary,  particles  can  be  grouped  into  several 
categories  based  not  strictly  on  size  but  on  the  value  of  ^  discussed  in  chapter  2.  Particles 
for  which  P  >  0.5  are  unbound  and  quickly  leave  the  system  on  hyperbolic  orbits.  These 
particles  are  referred  to  as  meteoroids.  Particles  that  are  below  but  still  close  to  the  blow 
out  limit  of  /3  =  0.5  are  called  ^  critical  particles  and  have  values  of  that  fall  in  the  range 
of  0.1  <  P  <  0.5.  These  particles  are  still  bound  but  have  their  orbits  significantly  affected 
by  radiation  effects.  In  general  the  distribution  of  these  particles  will  extend  outwards  from 
that  of  the  bodies  that  that  created  them.  Particles  for  which  ^  is  very  small  are  generally 
the  larger  particles  that  are  relatively  unaffected  by  raxliation  affects. 

Another  class  of  particles  is  the  one  in  which  PR  drag  plays  a  significant  role.  These 
particles  are  the  ones  that,  because  of  PR  drag,  spiral  in  towards  the  star  before  being 
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Figure  7-4:  The  effective  optical  depth  assumed  for  Beta  Pictoris  (the  solid  line)  and  the 
threshold  value  for  which  evolution  due  to  PR  drag  becomes  insignificant  (the  dotted  dashed 
line). 


destroyed  by  collisions  or  evaporation.  The  value  of  /3  for  these  particles  is  in  the  range  of 
Ppr  <  <  0.1,  where  Ppr  is  defined  as  the  value  of  ^  for  which  the  collisional  lifetime  is 

the  same  as  the  PR  drag  lifetime  (see  chapter  2).  This  value  is  expressed  as 

Ppr  =  5000  Te//(r)\/(M0/M*)(r’/o©).  (7.1) 

If  the  disk  has  a  high  enough  optical  depth,  these  particles  will  be  destroyed  by  collisions 
before  they  experience  significant  inward  migration  which  will  result  in  this  category  not 
existing  in  such  disks  (Wyatt  et  al.,  1999).  The  threshold  value  for  the  effective  optical 
depth  is  determined  by  finding  where  Ppr  >0.5  and  is  given  by 

re//(r)  >  10-V(^*/^©)(«®A)-  (7-2) 


Artymowicz  and  Clampin  (1997)  proposed  an  analytic  form  of  the  normal  optical  depth 
in  Beta  Pictoris  based  on  the  observed  surface  brightness.  This  expression  has  the  form 

1-1 


r(r)  =  2t„ 


\  —Pinner  /  ^ 

^m/  \^m 


^  \  pouter 


(7.3) 
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where  is  the  maximum  normal  optical  depth  which  is  reached  at  a  distance  of  r^.  The 
exponents  Pinner  and  Pouter  represent  the  slopes  of  the  inner  and  outer  disc.  The  values 
of  the  parameters,  as  estimated  from  12  pm  observations,  are  as  follows:  Tm  =  8.5  •  10“®, 
I'm  ~  50  AU,  Pinner  ~  2  and  Pouter  ~ 

Plots  of  the  assumed  optical  depth  and  the  threshold  value  are  shown  in  figure  7-4. 
This  figure  suggests  that  the  scenario  described  above  does  apply  and  that  the  disk  of  Beta 
Pictoris  does  not  possess  a  significant  distribution  of  PR  drag  affected  particles.  In  this 
situation,  half  of  the  cross  sectional  area  is  thought  to  come  from  the  /3  critical  particles 
(Wyatt  et  ah,  1999).  Though  the  value  of  ^  is  dependent  on  the  material  properties  of  the 
dust  (see  figure  2-1  in  chapter  2),  for  the  materials  considered  likely  to  constitute  the  dust 
in  Beta  Pictoris,  the  lower  limit  for  the  ^  critical  particles  corresponds  to  £>  ss  lO/xm. 

Typical  particle  sizes  determined  from  past  observations  and  analysis  were  thought  to 
have  a  radius  of  «  2  to  3  pm  (Artymowicz,  1997).  However,  this  estimate  was  done  prior  to 
the  release  of  the  Hipparcos  data,  which  revealed  a  greater  distance  to  Beta  Pictoris  than 
had  been  previously  thought  and,  consequently,  a  more  luminous  star.  Another  approach 
used  to  study  the  size  distribution  of  dust  within  the  Beta  Pictoris  system  was  a  theoretical 
study  using  a  proposed  collisional  model  of  dust  generation  (Krivov  et  ah,  2000).  This 
analysis  resulted  in  an  estimate  of  the  dominant  particle  size  to  be  ~  10  pm.  This  was 
based  on  the  result  that  showed  the  largest  contributor  to  cross  sectional  area  was  in  the 
range  just  below  the  blow  out  limit. 

For  the  modeling  presented  in  this  research,  the  typical  particle  size  was  assumed 
to  take  the  value  of  ~10  pm  in  diameter  for  dust  grains  composed  of  organic  refractory 
material.  This  corresponds  to  particles  slightly  larger  than  the  blow-out  limit  of  ^0  ^  0.5. 
The  justification  for  using  the  specific  material  composition  is  discussed  in  the  next  chapter. 

As  a  final  note,  it  is  important  to  keep  in  mind  that  though  radiation  effects  play  a 
significant  role  in  the  radial  distribution  of  dust,  the  orbital  planes  are  unaffected.  This  is 
important  because  this  research  focuses  on  warping  which  is  directly  related  to  the  mean 
orbital  planes  of  the  dust  particles  and  source  bodies.  Therefore,  even  if  the  estimate  on 
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particle  size  is  not  completely  accurate,  only  the  radial  distributions  are  affected  and  anal¬ 
ysis  of  the  symmetry  planes  is  still  arguably  valid. 

7.3.2  Timescales  of  Physical  Processes  in  Beta  Pictoris 

The  expressions  for  timescales  associated  with  each  of  the  primary  physical  processes 
within  the  disk  discussed  in  chapter  2  were  applied  to  Beta  Pictoris.  Using  the  parameters 
taken  from  observations,  the  timescales  were  determined  and  are  plotted  as  a  function 
of  semi-major  axis  in  figure  7-5.  These  calculations  were  most  dependent  on  physical 
characteristics  of  the  disk  such  as  mass/luminosity  of  the  central  star,  typical  particle  size, 
and  optical  depth. 

The  mass  and  luminosity  of  the  central  star  were  taken  from  the  results  based  on  the 
Hipparcos  data  (Crifo  et  al.  1997):  M  =  l.SM©  and  L  =  8.7Lq.  As  pointed  out  in  chapter 
5,  the  most  recent  estimate  determined  by  observing  the  stars  constituting  the  Beta  Pictoris 
moving  group  suggests  an  age  of  ~  20  Myr  (Navascues  et  al.  1999). 

Based  on  the  values  of  optical  depth  given  earlier,  the  collisional  timescales  are  much 
shorter  than  the  PR  lifetimes.  These  calculations  are  based  on  dust  modeled  as  particles 
composed  of  organic  refractory  material,  and  typical  particle  sizes  ranging  from  1.0  to  25.0 
microns.  The  associated  timescales  axe  plotted  in  figure  7-5. 

Analysis  of  the  different  timescales  suggests  once  again  that  the  collisional  lifetime  is 
the  shortest  of  the  timescales  and  thus  collisional  processes  dominate  over  all  others.  This 
implies  that  the  observed  dust  is  arguably  a  direct  tracer  of  the  larger,  secularly  perturbed, 
bodies.  This  supports  the  idea  that  by  modeling  the  distribution  of  bodies  affected  by 
secular  pertiurbations,  the  distribution  of  dust  within  the  system  can  be  reasonably  simu¬ 
lated.  However,  two  questions  must  first  be  addressed:  what  is  the  magnitude  of  the  secular 
timescale  and  how  does  it  compare  to  the  estimated  age  of  the  disk? 

The  appropriate  timescale  associated  with  the  secular  perturbations  is  the  proper  pro¬ 
cessional  rate  that  is  determined  from  calculating  both  the  A  and  B  values  of  equations 
3.7a  and  3.7c.  These  rates  can  be  converted  to  timescales  by  using  the  expression 


(7.4) 
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Figure  7-5:  The  timescales  associated  with  the  various  physical  processes  expected  within 
the  Beta  Pictoris  disk.  The  dotted  line  represents  the  Keplerian  period.  The  solid  line 
represents  the  collisional  timescales  for  the  particles  in  the  disk  that  make  up  most  of  the 
disk’s  cross  sectional  area.  The  three  dashed-dotted  lines  represent,  from  bottom  to  top, 
the  PR  timescales  associated  with  particles  with  D=l,  5,  25  microns  respectively. 


semi— major  axis  (AU)  semi— major  axis  (AU) 


Figure  7-6:  Proper  precessional  rate  for  a  four  planet  system  (left)  and  the  associated 
timescales  (right).  The  solid  hne  represents  the  results  for  one  of  the  systems  of  planets 
used  in  the  models  discussed  in  the  next  chapter.  The  dotted  line  represents  the  same 
system  but  with  the  mass  of  each  planet  scaled  by  0.5.  The  dotted  dashed  hne  in  the  left 
panel  represents  the  age  of  Beta  Pictoris.  The  solid  lines  represents  the  Keplerian  period. 
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Since  ^  is  a  function  of  semi-major  axis,  t  is  as  well.  A  plot  of  the  variation  of  secular 
timescales  is  shown  in  figure  7-6  where  the  A  values  are  associated  with  one  of  the  planetary 
configurations  used  in  the  modeling.  Also  shown  on  the  plot  for  reference  is  the  Keplerian 
timescale  and  the  estimated  age  of  the  Beta  Pictoris  system.  Plotted  in  these  diagrams 
are  the  results  firom  the  five  planet  planetary  configuration  used  in  the  models  presented 
in  the  next  chapter  and  the  results  for  the  same  system  using  a  mass  scale  factor  of  0.5  for 
each  of  the  planets.  This  figure  provides  insight  into  the  magnitude  of  change  that  would 
occur  in  the  secular  timescales  for  a  far  less  massive  system.  The  important  characteristic 
of  this  illustration  is  that  there  exists  plausible  configurations  in  which  secular  timescales 
are  in  fact  shorter  than  the  estimated  age  of  the  disk.  Furthermore,  this  statement  is  fairly 
robust  as  the  results  are  not  strongly  affected  by  significant  changes  in  the  total  mass  of  the 
planets  comprising  the  system.  This  answers  the  question  of  whether  or  not  the  disk  has  ex¬ 
isted  long  enough  to  evolve  under  the  influences  of  the  secular  perturbations.  It  should  also 
be  noted  that  the  20  Myr  age  is  an  estimate  of  the  minimum  age  of  the  Beta  Pictoris  system. 

7.3.3  Dust  Grain  Composition 

Prior  to  1998,  it  was  generally  accepted  that  the  dust  grains  within  the  Beta  Pictoris 
disk  belonged  to  one  of  two  populations:  within  ~  100  AU  the  material  was  considered 
to  likely  be  silicates  while  outside  ~  100  AU  the  material  was  considered  to  likely  be  ice 
(Artymowicz  1997).  Li  and  Greenberg  (1998)  studied  the  material  properties  in  more  detail 
and  suggested  that  the  dust  particles  in  the  inner  region  are  fluffy  aggregates  of  sihcate 
core/organic  refractory  mantle  grains.  In  the  outer  regions  these  aggregates  also  possess 
an  ice  mantle.  A  transition  in  material  properties  or  distribution  is  suggested  in  the  UF 
observations  by  the  change  in  slope  of  the  brightness  profile  at  ~  80  -  100  AU,  as  shown  in 
figure  7-7. 

The  primary  goal  of  the  simulations  performed  in  this  research  were  to  address  the 
effects  of  gravitational  perturbations  on  the  disk  morphology,  i.e.  the  spatial  distribution  of 
the  dust.  Therefore,  a  simplistic  approach  was  taken  in  selecting  the  material  composition 
of  the  dust  particles  in  the  model.  This  was  based  on  determining  which  of  the  dust  model 
candidates  produced  approximately  the  same  ratios  between  the  10  and  18  /i/m  emissions 
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Figure  7—7:  The  radial  brightness  profiles  along  the  midplane  in  a  5  AU  swath.  This  is 
plotted  on  a  log-log  scale  to  facilitate  a  linear  fit  and  determine  the  index  for  an  assinned 
powerlaw.  Three  regions  were  identified  for  the  line  fit.  The  boundaries  of  the  regions  and 
associated  power  law  indices  are  shown  on  the  plot. 
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as  determined  from  the  observations.  Implied  in  this  assumption  is  that  the  same  particles 
are  responsible  for  the  10  and  18  //m  emissions,  over  the  entire  disk. 

In  a  collisional  model,  the  grains  are  likely  to  be  porous  aggregates  of  silicate  core  - 
organic  refractory  mantle  particles  with  porosity  on  the  order  of  ^0.5  (Thebault  et  ah,  2003) 
as  opposed  to  cometary  dust  that  could  have  porosities  ~0.95  (Li  and  Greenberg,  1998).  For 
comparison,  data  from  dust  attributed  to  the  collisional  destruction  of  asteroids  in  the  solar 
system  suggests  a  porosity  of  0.50.  If  this  research  focused  on  an  exact  reproduction  of 
the  brightness  distribution,  a  more  thorough  analysis  would  be  required.  However,  since  the 
morphology  of  the  disk,  mostly  at  a  single  waveband,  is  the  primary  concern  of  the  modeling, 
the  analysis  considered  only  homogeneous  Mie  spheres  of  astronomical  silicates  and  organic 
refractories.  Though  aggregates  of  silicate  core  -  organic  refractory  mantles  are  likely 
present,  they  are  much  more  difficult  to  handle  through  the  simple  approach  used  in  this 
research.  In  order  to  handle  these  composite  materials  the  use  of  an  effective  medium  theory, 
such  as  Maxwell  Garnett,  is  necessary.  However,  Several  important  limitations  are  inherent 
to  this  approach.  First,  Maxwell-Garnett  effective  medium  theory  has  many  limitations  that 
are  discussed  thoroughly  in  Gustafson  et  al  (2001).  Among  these  limitations  is  the  necessity 
of  the  particles  making  up  the  aggregate  to  be  small  compared  to  the  wavelength  of  interest. 
Though  the  particles  making  up  an  aggregate  of  10  /im  with  a  fifty  percent  porosity  may 
fit  this  criteria  for  11  and  18  ^m  observations,  it  still  causes  concerns  when  considering 
the  fact  that  the  interaction  with  the  shorter  wavelengths  is  primarily  responsible  for  the 
heating  of  the  particles.  Furthermore,  Gustafson  et  al.  showed  significant  weaknesses  when 
considering  regions  near  the  10  //m  silicate  feature. 

Since  a  detailed  investigation  into  the  potential  material  composition  and  structure 
of  the  dust  grains  within  the  disk  was  well  outside  the  scope  of  this  research,  the  simple 
approach  presented  here  was  deemed  satisfactory.  To  reiterate,  if  a  detailed  analysis  of  the 
brightness  distribution  was  the  focus  of  the  research  and  not  disk  morphology  at  a  single 
waveband,  a  much  more  rigorous  and  thorough  analysis  would  be  required. 
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Temperature  profile  of  dust  particles 

To  simulate  the  emissions  from  dust  particles  it  is  necessary  to  make  an  estimate  of  the 
temperature  of  dust  grains  as  a  function  of  distance  from  the  star.  The  temperature  profiles 
were  determined  by  considering  the  scattering  efficiencies,  as  determined  from  Mie  theory, 
for  the  materials  thought  to  most  likely  represent  the  dust  grains  within  the  disk.  The 
complex  refractive  indices  and  the  temperature  profiles  for  the  different  materials  considered 
are  shown  in  figure  7-8.  This  diagram  lists  the  optical  constants  as  a  function  of  wavelength 
and  the  temperature  profiles  as  a  function  of  distance  from  the  star.  The  values  for  the 
optical  constants  for  the  astronomical  silicates  and  organic  refractories  were  taken  from 
Laor  and  Draine  (1993)  and  Li  and  Greenberg  (1997)  respectively.  The  results  for  several 
particle  diameters  are  shown. 

To  select  which  profile  is  appropriate  to  use  in  the  models,  the  thermal  emissions  at  10 
and  18  nm  were  compared  to  the  observations.  Figure  7-9  shows  the  temperature  and  flux 
ratios  as  a  function  of  size  for  various  radial  distances  from  the  star.  The  observed  mean 
flux  ratio  is  shown  as  a  horizontal  line  and  the  shaded  region  shows  the  one  sigma  uncer¬ 
tainty.  Only  the  organic  refractory  material  produces  reasonable  values  that  are  consistent 
with  the  previously  assumed  particle  size  of  ~10 


7.3.4  Mass  of  the  Disk 


The  total  mass  of  dust  within  a  disk  can  be  estimated  from  the  millimeter  flux  that  is 
primarily  associated  with  the  disk’s  gas.  The  fact  that  the  dust  material  is  not  detected  in 
these  observations  is  generally  accounted  for  through  the  use  of  an  estimated  gas  to  dust 
ratio.  Beckwith  (1999)  gave  an  expression  for  disk  mass 


=  0.03  Me  (j^)'  (^)' 


50  K  0.02  cm^g 


2„-l 


<  T  > 


^1.3mm 


(7.5) 


where  is  the  measured  flux  density,  D  is  the  distance  in  parsecs,  A  is  the  observational 
wavelength,  <  T  >  is  the  estimated  mean  temperature,  and  k  is  the  particle  mass  opacity 
which  is  ~  0.02(1.3  mm/A)cm^  G“^  (Beckwith,  1999).  The  millimeter  flux  density  for  the 
Beta  Pictoris  disk  was  measured  by  Liseau  et  al.  (2003)  at  24.3  ±  0.3  mJy  at  1200  nm. 
This  provides  an  estimated  mass  for  the  disk  as  ~  IM®.  This  is  consistent  with  the  values 
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wavelength  (/xm)  Distance  (AU) 

Figure  7-8:  Complex  index  of  refraction  (left  column)  for  several  materials  considered  likely 
candidates  for  the  dust  in  the  Beta  Pictoris  disk.  Also  shown,  is  the  temperature  profile 
(right  column)  as  a  function  of  distance  from  the  central  star. 
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Figure  7-9:  Temperataure  profile  (left  column)  for  dust  particles  as  a  function  of  particle 
size  in  the  Beta  Pictoris  system.  Also,  the  ratio  of  10  and  18/xm  flux  densities  (right  column) 
as  a  function  of  particle  size.  The  materials  are  the  same  as  those  considered  in  flgure  7-8. 
The  horizontal  line  in  the  plots  of  flux  ratio  show  the  observed  mean  value  and  shaded 
region  shows  the  one  sigma  uncertainty. 
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of  0.3  to  0.5  Mq  for  the  total  mass  of  dust  determined  by  Li  and  Greenberg  (1998)  and 
Artymowicz  (1997). 

The  mass  of  the  disk  is  important  to  consider  because  of  the  potential  influence  of  the 
disk  on  the  other  massive  bodies  and  dust.  As  will  be  discussed  in  the  following  sections, 
the  proposed  planetary  system  consists  of  planets  that  have  masses  comparable  to  that  of 
Jupiter,  ~  317  M®.  This  suggests  that  the  mass  of  the  disk,  which  is  spread  out  over  an 
enormous  area,  should  have  little  effect  on  the  dynamics  when  compared  to  that  of  the 
proposed  planets.  However,  it  needs  to  be  pointed  out  that  the  mass  associated  with  the 
source  bodies  is  not  considered  in  calculations  based  on  this  approach.  This  could  lead  to 
a  significant  underestimate  of  the  total  mass.  A  further  analysis  of  the  potential  influence 
of  a  massive  disk  is  presented  in  chapter  9. 


7.3.5  Distribution  of  Dust 

As  described  in  the  previous  chapter,  the  modeling  process  requires  an  estimate  of  the 
radial  distribution  of  orbits.  This  is  obtained  by  using  the  estimated  number  density  as 
estimated  from  observational  data.  As  discussed  earlier,  it  is  assumed  that  the  particles 
primarily  responsible  for  the  observed  thermal  emissions  are  rather  quickly  destroyed  by 
collisions  before  being  removed  from  the  system  by  radiation  pressure.  This  implies  that 
the  particles  not  only  need  a  source  for  replenishment,  but  that  the  observed  particles 
trace  the  distribution  of  source  bodies.  Prom  analysis  of  optical  depth  as  determined  from 
observations  and  from  theoretical  analysis  based  on  a  collisional  model  of  Beta  Pictoris,  an 
expression  for  the  number  density  was  given  by  Krivova  et  al.  (2000)  and  has  the  form 


where  the  values  of  pinner,  P&ater-,  and  were  suggested  to  be  roughly  1.0,  2.7,  and  60  AU 
respectively. 

The  important  information  provide  in  this  expression  is  the  general  shape  of  the  distri¬ 
bution  function.  The  models  in  this  research  require,  as  an  input  parameter,  a  normalized 
distribution  function  which  means  that  the  magnitude  of  no  is  unimportant.  The  actual 
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values  of  density  are  manifested  in  the  modeling  process  by  specifying  the  total  disk  surface 
area,  which  is  a  separate  model  parameter. 

7.3.6  The  Proper  Elements 

Necessary  for  the  simulation  of  a  disk  are  the  proper  orbital  elements;  these  parameters 
link  the  model  to  the  initial  conditions  associated  with  the  dust  at  the  moment  of  creation. 
The  values  used  for  the  models  in  this  research  are  based  on  the  unbiased  distribution 
determined  for  the  main-belt  asteroids  (Xu  et  al.,  1993).  This  distribution  however,  was 
scaled  in  such  a  way  as  to  provide  mean  values  consistent  with  observable  quantities  of 
the  disk.  As  discussed  in  chapter  4,  the  half-opening  angle  of  the  disk  is  related  to  the 
proper  inclination.  The  observed  opening  angle  has  been  measured  to  be  ~12  degrees 
(Artymowicz  et  al.,  1989)  which  implies  a  mean  proper  inclination  of  ~6  degrees.  Because 
of  the  edge-on  view  of  the  Beta  Pictoris  disk,  it  is  more  difficult  to  make  estimates  of  the 
proper  eccentricity.  However,  a  value  for  the  mean  proper  eccentricity  of  0.1  was  chosen 
in  order  to  maintain  the  same  relationship  between  the  mean  proper  inclination  and  mean 
proper  eccentricity  as  is  observed  for  the  main  belt  asteroids,  i.e.  <  sin/  >  «  <  e  >. 

The  probability  distribution  function  (pdf)  for  the  main  belt  asteroids  and  the  modified 
pdf  for  the  Beta  Pictoris  disk  are  shown  in  figure  7-10.  Also  shown  is  the  histogram 
representing  the  randomly  selected  values  of  proper  elements  based  on  the  modified  pdf. 

These  values  were  appropriate  based  on  the  values  for  the  half-opening  angle  docu¬ 
mented  in  the  literature.  However,  the  next  chapter  will  show  that  the  disk  may  actually 
be  better  represented  by  smaller  mean  values.  The  high  resolution  images  reveal  fairly  well 
isolated  regions  of  dust  that  are  suggestive  of  a  mean  proper  inclination  of  ~  3  degrees  and 
a  mean  proper  eccentricity  of  ~  0.05.  Several  models  were  run  with  these  values  and  did  in 
fact  produce  a  better  visual  match  to  the  features  of  the  high  resolution  images.  This  will 
be  discussed  further  in  the  following  chapter. 
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Figure  7-10:  Distribution  of  proper  orbital  elements  used  in  the  simulated  disk  of  Beta 
Pictoris.  The  solid  line  shows  the  theoretical  pdf  curve  determined  by  modifying  the  pdf  for 
the  main  belt  asteroids.  The  dotted-dashed  line  shows  the  pdf  for  the  main-belt  asteroids  as 
determined  from  the  distribution  of  Xu  et  al.  (1993).  The  histogram  shows  the  distribution 
of  values  randomly  drawn  using  the  modified  pdf  (solid  curve) . 

7.4  Summary 

This  chapter  provided  the  current  understanding  of  the  environment  around  Beta  Pic¬ 
toris  to  the  extent  necessary  to  create  adequate  models  for  an  investigation  into  disk  mor¬ 
phology.  The  next  chapter  will  discuss  the  way  in  which  these  characteristics  are  imple¬ 
mented  into  the  modeling  of  the  morphology  of  the  dust  disk  surrounding  Beta  Pictoris. 


CHAPTER  8 

APPLICATIONS  OF  SECULAR  PERTURBATION 
THEORY  TO  BETA  PICTORIS 

8.1  Introduction 

The  previous  chapter  discussed  the  characteristics  of  the  physical  environment  sur¬ 
rounding  Beta  Pictoris.  These  characteristics  provide  the  boundary  conditions  needed  to 
develop  a  reasonable  model  of  the  circumstellar  disk.  However,  the  important  characteris¬ 
tics  needed  to  assess  the  warping  of  the  midplane  of  the  disk  are  associated  with  the  details 
of  the  orbits  of  the  perturbing  planets.  Therefore,  an  approach  must  be  defined  to  extract 
information  from  the  observations  of  a  disk  that  can  be  used  to  characterize  the  orbits  of 
an  embedded  planetary  system. 

A  system  of  multiple  planets  will  be  utilized  to  explain  the  existence  and  configuration 
of  the  major  symmetry  planes  of  the  disk.  However,  current  theory  of  planet-disk  interaction 
suggests  that  these  planets  will  also  have  secondary  effects  on  the  disk,  some  of  which 
should  be  observable  in  the  observations.  These  effects  include  gap  clearing  and  warping 
on  smaller  scales  -  both  of  which  are  suggested  by  the  UF  observations.  Regions  will  be 
depleted  through  interaction  with  planets  via  tidal  forces  or  accretion.  Both  process  can  be 
related  to  characteristics  of  the  planetary  orbits  or  the  masses  of  the  planets  themselves. 
Furthermore,  the  eccentricity  imposed  on  the  dust  particles  through  secular  perturbations 
combined  with  the  particles’  inherent  eccentricity  can  cause  the  dust,  not  initially  in  a 
region  cleared  by  a  planet,  to  enter  an  orbit  that  will  cause  it  to  enter  such  a  region,  and 
thus  be  removed  from  the  disk.  This  potentially  results  in  additional  regions  in  which  dust 
depletion  should  occur  and  has  implications  for  the  stability  of  a  large  disk  in  a  multi-planet 
system.  The  small  scale  warping  is  also  expected  for  planets  with  mutual  inclinations  due 
to  the  potentially  dramatic  behavior  of  the  forced  orbital  elements,  particularly  the  forced 
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inclination  and  node.  By  probing  the  small  scales,  information  about  the  behavior  of  the 
forced  parameters  can  be  gleaned. 

These  quantifiable  features  provide  additional  constraints  on  the  perturbations  from 
the  planets  and  thus  on  the  planets  themselves.  The  way  in  which  these  constraints  are 
determined  and  then  apphed  is  the  primary  goal  of  this  chapter.  Sample  models  based  on 
the  approach  proposed  here  will  be  presented  along  with  a  discussion  of  the  results  in  the 
context  of  the  observations. 

8.2  Constraints  on  a  Planetary  System  from  Observations 

As  discussed  in  the  previous  chapter,  secular  perturbation  theory  predicts  that  the 
gravitational  influence  of  a  system  consisting  of  multiple  planets  will  define  the  mean  or¬ 
bital  planes  of  the  disk’s  constituent  dust  particles.  If  the  planets  have  a  nonzero  mutual 
inclination,  the  components  of  the  dust  particles’  orbital  elements  that  arise  from  the  grav¬ 
itational  perturbations,  the  forced  elements,  will  vary  with  semi-major  axis.  Since  the 
forced  inclination  and  longitude  of  ascending  node  define  the  planes  about  which  particles 
will  precess,  and  since  these  elements  are  a  function  of  semi-major  axis,  warping  occurs. 
This  effect  can  also  be  considered  in  terms  of  angular  momentum  planes.  At  radial  dis¬ 
tances  corresponding  to  the  semi-major  axis  of  one  of  the  planets,  the  dust  particles  will 
precess  around  the  same  plane  as  the  corresponding  planet  -  the  plane  perpendicular  to  the 
planet’s  angular  momentum  vector.  This  will  be  the  case  for  each  of  the  planets.  Between 
the  planets,  the  effect,  and  hence  the  behavior,  is  more  complicated  as  the  particles  are  now 
influenced  by  multiple  masses. 

The  choice  of  orbital  elements  of  a  proposed  system  requires  adherence  to  two  con¬ 
straints.  First,  the  chosen  planetary  configmration  must  have  an  associated  forced  inclina¬ 
tion  that  would  produce  the  large  scale  warping  well  documented  in  the  literature.  Second, 
the  secondary  effects  associated  with  the  planets  also  need  to  be  consistent  with  the  sections 
of  disk  identified  in  the  UF  images.  The  observations  of  Wahhaj  et  al.  (2003),  Weinberger 
et  al.  (2003),  and  Telesco  et  al.  (2003)  were  used  to  quantify  the  angle  of  the  inner  com¬ 
ponent  that  peaked  in  amplitude  at  ~20  AU.  The  observations  and  analysis  of  Heap  et  al. 
(2000)  were  used  to  determine  the  angle  between  the  intermediate  component  and  the  outer 
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component.  Since  the  UF  observations  were  the  highest  in  quality,  they  alone  were  used  to 
study  the  small  scale  features. 

The  following  sections  describe  how  the  details  of  the  observations  are  utilized  to  con¬ 
strain  the  characteristics  of  the  proposed  planetary  system.  First,  the  major  symmetry 
planes  are  used  to  constrain  the  major  angular  momentum  planes  within  the  disk.  Second, 
the  additional  features  observed  in  the  high  resolution  images  are  used  to  further  constrain 
the  models  by  suggesting  the  likely  locations  of  the  planets. 

8.2.1  Large  Scale  Symmetry  Planes 

Figures  8-1  and  8-2  compare  the  UF  observations  with  the  previously  identified  major 
symmetry  planes.  Figure  8-1  shows  a  contour  plot  of  both  the  18  and  11  pm  images.  The 
dashed  and  dotted  lines  indicate  the  major  symmetry  planes  of  the  inner  disk  (Weinberger 
et  al.,  2003)  and  the  larger  scale  ridge-line  (Heap  et  al.,  2000).  Figure  8-2  shows  the  same 
data,  but  each  vertical  slice  has  been  fit  by  a  gaussian.  Contours  are  shown  for  the  10% 
values  and  the  FWHM. 

The  importance  of  these  two  plots  is  the  consistency  between  the  previously  identified 
symmetry  planes  and  the  UF  data.  However,  though  the  images  are  consistent  with  the 
previous  assertions,  they  don’t  match  exactly.  The  inner  disk  is  only  tilted  by  8-9  degrees  as 
opposed  to  the  previous  estimate  of  12  degrees.  As  discussed  in  chapter  4,  major  symmetry 
planes  are  associated  with  primary  angular  momentum  planes  within  the  system.  The  inner 
symmetry  planes  are  associated  with  the  angular  momentum  planes  of  the  planets  and  the 
outer  large  scale  symmetry  plane  is  associated  with  the  total  system  angular  momentum 
plane.  These  facts  provide  constraints  on  the  masses  and  locations  of  any  planets  within 
the  system. 

These  constraints  can  be  assessed  by  considering  the  relationships  between  the  masses 
and  semi-major  axes  of  the  planets.  Figure  8-3  shows  the  results  from  a  calculation  of  the 
inclination  of  the  plane  associated  with  the  system’s  total  angular  momentum  with  respect 
to  the  plane  of  the  second  planet.  This  calculation  was  made  by  assuming  an  expression 
for  the  inclination  that  was  simplified  to  a  function  of  the  mass  ratios  and  the  ratio  of 
semi-major  axes.  To  make  this  simplification,  several  assumptions  were  made.  First,  the 
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Figure  8-2;  Same  as  figure  8-1  but  with  each  vertical  slice  normalized  to  the  peak  value 
and  fitted  by  a  gaussian.  Plotted  are  the  contoiurs  representing,  for  each  vertical  cut,  the 
>  99%  brightness,  the  FWHM,  and  the  10%  values.  These  are  also  shown  in  comparison 
to  the  STIS  and  Weinberger  reference  planes. 
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M2y/M  1 

Figure  8-3;  Calculated  values  of  the  inclination  of  the  plane  defined  by  the  total  angular 
momentum  of  the  system  as  a  function  of  the  mass  ratios  of  the  two  planets  assumed  to 
comprise  the  system.  Each  curve  represents  a  unique  value  of  the  semi-major  axis  ratio. 
The  observed  value  for  the  inclination  of  the  invariant  plane  is  shown.  The  intersection  of 
this  line  with  the  calculated  curves  defines  the  acceptable  values  of  semi-major  axis  and 
mass  ratios.  The  reference  plane  is  that  of  the  outer  planet. 
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Figure  8-4:  The  forced  inclination  for  two  planets  with  mass  and  semi-major  axis  ratios 
determined  from  the  plot  shown  in  figure  8-3.  The  filled  circles  indicate  the  inclinations 
of  the  planets.  The  reference  plane  is  that  of  the  outer  planet.  The  figure  on  the  left 
corresponds  to  the  02/01  =  2.5  and  the  figure  on  the  right  corresponds  to  case  of  02/01  =  3.5. 


planes  were  assumed  to  be  a  result  of  just  two  planets  -  one  associated  with  each  of  the 
inner  symmetry  planes.  Furthermore,  the  inclinations  of  each  planet  were  assumed  to  be 
the  same  as  the  inclinations  of  the  corresponding  symmetry  planes.  The  final  assumption 
was  that  the  nodes,  including  that  of  the  total  angular  momentum  plane,  were  aligned  with 
the  line  of  sight  of  the  observation. 

The  curves  shown  in  figure  8—3  shows  the  inclination  of  the  system’s  total  angular 
momentum  plane  as  a  function  of  mass  ratios  of  the  two  planets.  Each  curve  represents  a 
different  value  for  the  ratio  of  semi-major  axes.  Also  plotted  is  the  value  for  the  observed 
inclination,  ~  4.8  degrees.  It  is  the  intersection  of  this  line  with  the  various  curves  that 
provide  the  acceptable  value  for  the  mass  ratio  for  a  given  value  of  semi-major  axis  ratio. 

The  information  plotted  in  figure  8-3  can  be  used  to  create  a  model  system  with  major 
symmetry  planes  matching  those  of  Beta  Pictoris.  Examples  of  the  forced  inclinations 
associated  with  planetary  systems  defined  by  the  ratios  determined  from  figure  8-3  are 
shown  in  figure  8-4.  Also  plotted  in  this  figure  are  the  inclinations  associated  with  the 
planets  composing  the  system  -  the  filled  circles.  The  semi-major  axis  of  the  inner  planet 
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Figure  8-5:  A  model  disk  for  a  two-planet  system.  The  orbital  parameters  and  masses  used 
were  determined  from  the  calculation  shown  in  figure  8-3.  The  planetary  system  used  is 
the  same  system  that  produced  the  forced  inclination  shown  in  the  second  panel  of  figure 
8-4. 
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in  both  plots  is  15  AU  and  corresponds  to  the  estimated  radii  of  the  inner  hole.  The  semi¬ 
major  axis  ratio  is  2.5  in  the  first  plot  and  3.5  in  the  second  plot.  The  mass  ratios  in  the 
first  and  second  plots  are  0.47  and  0.40  respectively. 

Beyond  the  assumptions  already  discussed,  one  additional  constraint  needed  to  be 
considered  -  the  radial  extent  of  the  middle  or  intermediate  symmetry  plane.  Heap  et  al. 
(2000)  showed  the  central  disk  to  extend  out  to  about  70  AU  which  is  consistent  with  the 
value  estimated  from  the  UF  observations  of  about  70-80  AU.  This  feature  corresponds  to 
the  region  where  the  forced  inclination  passes  through  a  significant  secular  resonance  and 
begins  approaching  its  invariant  value.  This  can  be  seen  in  the  plots  of  figure  8-4.  For  the 
first  configuration,  this  occurs  at  approximately  50-60  AU.  For  the  second  configmation, 
this  occurs  at  approximately  80-90  AU.  To  show  how  this  translates  into  a  disk  model,  a 
simple  model  was  developed  using  the  second  of  the  configurations  shown  in  figure  8-4. 
A  visualization  of  this  simple  model  is  shown  in  figure  8-5.  In  this  diagram,  the  multiple 
components  can  be  clearly  seen  and  the  relationship  between  transitional  locations  and 
features  in  the  forced  inclination  are  obvious. 

It  should  be  noted  that  with  the  limited  information  available  this  is  as  far  as  an 
analysis  of  the  major  symmetry  planes  can  go.  Because  angulax  momentum  is  the  key 
issue,  it  is  only  the  ratios  of  mass  and  semi-major  axis  that  can  be  constrained  by  this 
approach.  However,  the  identification  of  additional  features  allows  a  deeper  probing  of  the 
system.  The  way  in  which  this  is  accomplished  is  discussed  in  the  next  section. 

8.2.2  Brightness  Enhancements  and  Voids 

A  different  reduction  of  the  UF  data  is  shown  in  figure  8-6.  This  image  shows  the 
rotated  18  /zm  image  with  irregular  contours  used  to  highlight  the  observed  structure. 
Along  the  mid-plane  of  the  disk  image,  multiple  brightness  enhancements  and  voids  were 
observed.  The  voids  are  arguably  logical  locations  for  planets  based  on  current  theories 
describing  the  accretion/clearing  process  of  planetary  formation  (Lin  and  Papaloizou  1993). 
As  was  discussed  in  chapter  2  the  region  cleared  during  the  process  of  planet  formation  is 
roughly  four  times  the  radius  of  the  Hill’s  sphere.  For  Beta  Pictoris,  M*  =  I.8M0,  this 
corresponds  to  regions  that  extend  approxiamately  r  «  0.224  a  (mp/Mj)^^^.  For  the  values 
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Figure  8-6:  (left)  18  /im  image  of  thermal  emissions  from  Beta  Pictoris. 


of  mass  and  semi-major  axis  considered  in  this  research,  this  corresponded  to  roughly  5-6 
AU.  These  values  are  consistent  with  the  typical  distances  between  the  bright  regions  which 
were  approximately  6-10  AU. 

Gaps  are  generally  attributed  to  tidal  interactions  between  the  planets  or  protoplanets 
and  the  disk.  As  a  result  of  these  interactions,  the  material  interior  to  the  planet’s  orbit 
drifts  inward  and  the  material  exterior  to  the  planet’s  orbit  drifts  outward.  The  inward 
drift  of  material  is  a  result  of  angular  momentum  being  transferred  from  the  disk  to  the 
planet  and  the  outward  drift  of  material  is  a  result  of  angular  momentum  being  transferred 
from  the  planet  to  the  disk.  Though  this  work  will  not  investigate  the  possibility,  resonant 
trapping  may  also  play  a  significant  role  in  creating  regions  that  appear  as  bright  spots  in  the 
observations.  As  is  the  case  for  the  asteroid  belt,  mean  motion  and/or  secular  resonances 
are  likely  to  play  roles  in  defining  the  boundaries  for  different  reservoirs  of  asteroid  like 
bodies. 

The  UF  images  of  Beta  Pictoris  indicate  regions  in  which  gaps  within  the  disk  may 
exist.  Brightness  enhancements  were  identified  in  the  observations  and  are  listed  in  table 
8-1.  The  regions  between  these  bright  spots  are  likely  depleted  regions  and  are  arguably 
related  to  a  gap  clearing  process.  Based  on  the  assumption  that  these  potential  gaps  are  a 
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Table  8-1:  Brightness  Peaks 


NE  Wing 

SW  Wing 

Radial 

Position  Relative 

Radial 

Position  Relative 

Position 

to  mid-plane 

Position 

to  mid-plane 

26  AU 

North 

26  AU 

South 

32  AU 

South 

36  AU 

North 

43  AU 

North 

42  AU 

South 

49  AU 

South 

52  AU 

North 

58  AU 

North 

58  AU 

South 

68  AU 

North 

68  AU 

South 

79  AU 

North 

77  AU 

South 

93  AU 

North 

96  AU 

South 

result  of  clearing  by  planetary  bodies,  the  corresponding  radii  were  used  as  possible  ranges 
for  the  semi-major  axes  of  the  planets  constituting  the  model  system. 

The  gaps  may  also  provide  a  reasonable  means  for  estimating  the  eccentricities  of  the 
orbits.  Based  on  the  analysis  done  by  Sari  and  Goldreich  (2003)  there  is  a  symbiotic  re¬ 
lationship  between  the  disk  and  the  massive  planets  within.  This  relationship  provides  a 
connection  between  the  eccentricity  of  the  planets  and  the  gaps  which  they  clear.  Sari  and 
Goldreich  suggest  that  once  the  gap  is  cleared,  the  eccentricity  of  the  planet  is  maintained 
at  a  value  of  e  «  w/r,  where  e  is  the  planet’s  eccentricity,  w  is  the  gap  width,  and  r  is  the 
distance  from  the  star.  Therefore,  the  values  from  the  18  fim  image,  where  the  gaps  are 
most  clearly  observed,  provide  a  constraint  to  the  final  model.  However,  several  problems 
are  associated  with  determining  the  values  of  eccentricity  with  this  approach.  First,  the 
exact  edges  of  the  inferred  gaps  are  difiicult  to  determine.  Second,  even  a  conservative 
estimate  of  the  edge  locations  suggest  extremely  large  eccentricities  for  the  planets,  some 
as  high  as  e  ~  0.3,  which  are  not  consistent  with  estimates  based  on  other  considerations 
such  as  those  presented  in  the  next  section.  Therefore,  the  eccentricities  determined  from 
this  type  of  analysis  of  the  brightness  voids  will  be  used  only  as  an  absolute  upper  limit. 

8.2.3  Forced  Orbital  Elements 

The  bright  spots  described  in  table  8-1  hint  at  a  symmetry  in  their  locations  from  the 
central  star  and  in  their  relative  position  to  the  mid-plane  of  the  mid-IR  images.  Though 
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not  conclusively,  this  suggests  a  relationship  between  regions  in  the  NW  extension  and  in 
the  SW  extension.  If  these  regions  are  related,  they  can  be  interpreted  as  toroidal  sections  of 
disk,  similar  to  those  described  in  chapter  4.  If  this  is  true,  they  should  have  characteristics 
related  to  the  forced  and  proper  orbital  elements  as  described  in  that  chapter  and  illustrated 
in  figure  4-6.  This  implies  that  the  offset  of  the  center  of  symmetry  firom  the  star  for  each 
section  of  disk  should  relate  to  the  local  value  of  the  forced  eccentricity.  Similarly,  the 
inclination,  with  respect  to  a  reference  plane,  should  relate  to  the  local  value  of  the  forced 
inclination  with  respect  to  the  same  reference  plane. 

Extracting  the  offset  and  inclination  from  the  observations,  for  a  given  value  of  semi¬ 
major  axis,  allows  the  development  of  a  tabulated  set  of  forced  elements.  Though  these  do 
not  strictly  define  the  form  of  the  forced  orbital  elements,  they  do  set  constraints  to  which 
a  final  form  of  the  orbital  elements  must  adhere. 

The  way  in  which  the  information  about  the  local  values  of  the  forced  elements  was 
extracted  was  by  reproducing  the  sections.  This  was  accomplished  by  reproducing  the 
different  sections  of  disk  by  defining  a  set  of  forced  and  proper  orbital  elements  for  each 
section.  The  proper  elements  were  chosen  in  a  way  consistent  with  the  relationships  dis¬ 
cussed  in  chapter  4:  the  radial  width  of  the  section  is  ~  2aep  and  the  angular  thickness 
is  ~  2ip.  Likewise,  the  forced  elements  were  chosen  in  a  way  consistent  with  the  offset 
from  the  center  of  symmetry  being  ~  ae^  and  the  inclination  of  the  section  of  disk  from  a 
reference  plane  being  ~  i/. 

Using  the  parameters  that  provided  the  best  visual  fit,  the  sections  of  disk  consid¬ 
ered  were  reproduced  and  are  shown  in  figure  8-7.  The  figure  shows  first  the  sections  of 
disk  with  no  inclination,  but  with  an  eccentricity  and  orientation  (longitude  of  pericenter) 
capable  of  reproducing  the  offsets  from  the  center  of  symmetry  as  observed  from  the  obser¬ 
vational  line-of-sight.  It  should  be  emphasized  that  these  two  parameters  are  coupled  and 
cannot  be  separated  with  only  a  single  line-of-sight.  Second,  the  same  sections  axe  shown 
but  with  an  inclination  and  orientation  (longitude  of  ascending  node)  incorporated  that  is 
consistent  with  the  inclination  from  the  reference  mid-plane  of  the  mid-IR  images.  The 
final  two  panels  show  the  edge-on  view  of  the  system  with  the  final  panel  only  showing  the 
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Figure  8-7;  Reproduction  of  the  sections  of  disk  possibly  identified  from  the  UF  observa¬ 
tions.  a)  disk  sections  with  no  inclination  but  with  an  eccentricity  capable  of  producing  the 
observed  offsets  in  symmetry,  b)  same  as  panel  a  but  with  a  mutual  inclination  capable  of 
producing  the  observed  positions  of  the  brightness  enhancements  relative  to  the  mid-plane, 
c)  a  side  view  of  the  configuration  shown  in  panel  b.  d)  the  same  as  panel  c  but  with  only 
the  regions  of  greatest  optical  depth  shown. 
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Figure  8-8:  A  plot  of  the  orbital  elements  associated  with  the  disk  sections  shown  in  figure 
8-7.  Also  shown  is  a  polynomial  fit  to  the  points  associated  with  each  disk  section.  Though 
this  does  not  define  the  functional  form  of  the  forced  elements  with  semi  major  axis,  it  does 
show  the  constraints  to  which  any  solution  must  conform.  The  grey  areas  correspond  to 
the  regions  in  which  planets  are  most  likely  to  be  found. 
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Figure  8-9:  Simulated  disk  based  on  a  set  of  forced  elements  extracted  from  the  UF  images. 
The  values  used  for  the  forced  elements  are  those  plotted  in  figure  8-8. 


regions  corresponding  to  those  regions  with  greatest  optical  depth  -  the  regions  that  would 
correspond  to  the  observed  bright  spots  in  the  observations. 

The  values  for  the  parameters  were  found  by  adjusting  the  eccentricities  and  pericenters 
until  the  distances  from  the  star  to  the  bright  regions,  in  both  extensions,  matched  those 
reported  in  table  8-1.  The  inclinations  and  nodes  were  adjusted  until  the  orientation,  with 
respect  to  the  mid-plane,  for  each  bright  spot  roughly  matched  those  seen  in  the  image 
shown  in  figure  8-6. 

The  values  of  the  forced  parameters  that  were  used  to  reproduce  the  locations  of  the 
bright  spots  were  tabulated  and  are  shown  graphically  in  figure  8-8.  This  figure  not  only 
shows  the  complete  set  of  forced  orbital  parameters,  but  it  also  shows  the  regions  in  which 
planets  are  most  likely  to  be  found  -  the  shaded  regions  in  the  diagram.  This  judgment 
is  based  on  the  simple  idea  that  a  planet  will  clear  a  region  of  the  disk,  and  thus  produce 
a  gap.  However,  it  should  be  noted  that  because  a  planet  will  create  a  gap,  it  cannot  be 
automatically  assumed  that  the  apparent  presence  of  a  gap  is  conclusive  evidence  of  a  planet. 
It  is  also  possible  that  the  brightness  enhancements  are  in  fact  caused  by  resonant  trapping 
from  the  gravitational  influences  from  the  planets.  Secular  resonances  may  also  play  a  role  in 
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Figure  8-10:  Superposition  of  the  estimated  brightness  locations  from  figure  8-8  over  the 
18  //m  data. 
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Table  8-2:  Planet  Locations  and  Masses 


Planet 

a  (AU) 

Mass  (10  M*.) 

Mftss 

1 

21.2 

11.20 

2.12 

2 

31.0 

2.50 

0.47 

3 

38.4 

1.25 

0.24 

4 

46.0 

0.50 

0.10 

5 

55.5 

0.30 

0.06 

6 

63.0 

0.25 

0.05 

7 

73.5 

0.25 

0.05 

defining  the  edge  of  of  a  distribution  of  source  bodies,  much  as  they  due  in  the  asteroid  belt 
of  the  solar  system.  However,  if  it  is  assumed  that  the  gaps  are  associated  with  planets,  the 
brightness  enhancements  can  be  logically  used  to  assess  the  local  values  of  the  forced  orbital 
elements.  To  ensure  that  a  set  of  orbital  elements  consistent  with  these  values  will  in  fact 
reproduce  the  observed  features,  a  model  disk  was  simulated  using  the  extracted  parameters 
as  the  forced  elements.  The  result  is  shown  in  figure  8-9.  This  image  reproduces  the  location 
of  the  brightness  enhancements  and  more  importantly  the  orientation  of  each  enhancement 
with  respect  to  the  mid-plane.  This  validates  the  assumptions  made  in  extracting  the 
estimates  of  local  values  for  the  forced  orbital  elements  from  the  observations.  Figure  8-10 
shows  a  comparison  of  the  locations  of  the  bright  spots  in  figure  8-7  to  the  18  /tfm  images 
of  Beta  Pictbris. 

Once  again,  it  is  important  to  point  out  that' the  inclination  and  node  are  likely  better 
constraints  than  the  eccentricity  and  pericenter.  Though  both  are  coupled,  node  to  inclina¬ 
tion  and  eccentricity  to  pericenter,  the  later  set  is  much  more  difficult  to  constrain  because 
of  the  single  line  of  sight.  This  fact  will  be  important  in  the  following  step  of  developing  a 
proposed  planetary  system. 

The  initial  configuration  of  planets  used  was  based  on  the  assmnption  that  each  region 
identified  as  a  prospective  gap  corresponded  to  a  plant,  and  thus  set  the  value  for  the  semi¬ 
major  axis.  Figure  8-11  shows  a  complete  set  of  forced  orbital  elements  calculated  from 
secular  perturbation  theory  as  described  in  chapter  3.  This  set  is  based  on  seven  planets, 
represented  by  the  filled  circles,  placed  within  the  regions  identified  as  likely  locations  for 
planets  (the  shaded  regions).  This  model  is  constrained  by  the  values  determined  from 
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Figure  8-11:  A  complete  set  of  forced  orbital  elements  as  a  function  of  semi-major  axis  for 
a  seven  planet  system.  The  masses  associated  with  the  planets  are  shown  in  table  8-2. 
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analysis  of  the  disk  sections  (shown  as  diamonds  on  the  plots)  and  by  the  large  scale 
symmetry  planes. 

The  approach  in  selecting  a  set  of  orbital  elements  took  two  steps.  The  first  step  was  to 
investigate  whether  or  not  a  system  of  planets  located  in  each  of  the  regions  interpreted  as 
gaps  could  produce  a  set  of  forced  orbital  elements  consistent  with  the  constraints  extracted 
firom  the  observations  -  the  observed  sections  of  disk.  The  secondary  goal  was  to  investigate 
whether  or  not  a  fewer  number  of  planets  could  also  produce  a  set  of  orbital  elements  that 
was  consistent  with  the  constraints  described  above.  The  absolute  minimum  number  of 
planets  needed  was  based  on  the  number  of  major  symmetry  planes.  Since  there  are  three 
well  documented  symmetry  planes,  a  minimum  of  two  planets  was  needed  (see  discussion 
in  chapter  4).  However,  it  was  found  through  trial  and  error  that  the  behavior  of  the 
forced  inclinations,  as  illustrated  in  figure  8-8,  could  not  be  reproduced  without  additional 
masses.  Planets  were  added  to  the  gap  locations,  one  at  a  time,  until  the  calculated  forced 
elements  provided  a  reasonable  match  to  the  constraints  extracted  from  the  observations. 
The  locations  of  the  planets  were  chosen  based  on  the  details  of  the  warp.  The  specific  gap 
locations  were  chosen  by  determining  where  the  most  dramatic  changes  occurred  between 
two  consecutive  regions  of  disk. 

The  Tninimnm  number  of  planets  is  clearly  difficult  to  determine.  The  number  is  clearly 
more  than  two,  but  beyond  that,  the  number  becomes  difficult  to  constrain.  An  example 
of  a  calculation  of  forced  orbital  parameters,  highlighting  the  problem  in  determining  the 
minimum  number  of  planets,  is  shown  in  figure  8-12.  This  figure  shows  the  forced  inclina¬ 
tion  and  node  for  a  two  and  three  planet  system.  Also  shown  on  the  plot  are  the  values 
extracted  from  the  observations.  The  two  planet  system  clearly  cannot  produce  the  com¬ 
plex  structure  suggested  by  the  observations.  The  three  plajiet  system,  on  the  other  hand, 
at  first  inspection  hints  at  the  possibility  of  fitting  the  observed  values  of  forced  inclination. 
However,  when  the  calculated  and  observed  values  for  the  forced  longitude  of  ascending 
node  are  examined,  the  three  planet  system  clearly  is  unable  to  reproduce  the  necessary 
values  of  the  forced  node.  This  is  significant  because  it  is  the  value  for  the  forced  node 
that  defines  the  orientation  of  the  observed  sections  of  disk.  Through  further  iterations  of 
placing  and  moving  planets,  it  was  found  that  a  five  planet  system  could  also  reproduce  the 
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Figure  8-12:  The  forced  inclination  and  node  for  a  two  (top)  and  three  (bottom)  planet 
system.  The  diamonds  correspond  to  the  values  extracted  from  the  observational  data.  The 
masses  for  the  planets  are  Mi  =  O.SMj  and  M2  =  0.2Mj  in  the  top  set  and  Mi  =  0.7Mj, 
M2  =  O.2M7,  and  Ms  =  0.2Mj  in  the  bottom  set. 
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Figure  8-13:  A  complete  set  of  forced  orbital  elements  as  a  function  of  semi-major  axis  for 
a  five  planet  system. 

necessary  forced  elements.  The  forced  orbital  elements  calculated  for  a  five  planet  system 
are  shown  in  figure  8-13.  This  provides  a  further  example  of  how  a  multiple  planet  system 
is  capable  of  meeting  the  constraints  set  by  analysis  of  the  observed  sections  of  disk.  This 
type  of  analysis  would  therefore  suggest  that  the  number  of  planets  is  between  five  and 
seven. 

However,  one  important  characteristic  about  the  form  of  the  forced  orbital  elements 
for  the  five  planet  system  must  be  made  clear.  Though  the  plot  of  the  forced  elements 
passes  through  the  points  corresponding  to  the  observed  sections  of  disk,  it  is  the  dramatic 
behavior  of  the  curve  that  allows  this  to  happen.  The  case  is  the  same  for  the  six  planet 
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systems.  What  this  implies  is  that  for  such  a  system  of  planets  to  produce  the  type  of 
morphology  observed  in  Beta  Pictoris,  the  high  amplitude  perturbations  are  required  to 
play  a  significant  role.  The  difficulty  with  this  is  that  a  region  aflfected  by  high  amplitude 
perturbations  is  likely  to  not  exist  for  a  long  period  of  time  and  thus  not  necessarily  relate 
to  truly  observable  features.  It  is  for  this  reason,  that  the  most  likely  candidate  is  the  seven 
planet  system. 

The  values  of  the  resulting  forced  elements  were  constrained  by  the  observed  disk 
sections  but  the  masses  were  primarily  constrained  by  angular  momentum  considerations. 
As  was  discussed  earlier,  the  large  scale  symmetry  planes  are  related  to  primary  angular 
momentum  planes  which  are  heavily  dependent  on  the  planetary  masses.  Ergo,  the  masses 
of  the  planets  were  adjusted,  after  the  locations  were  chosen,  until  the  large  scale  symmetry 
planes  matched  those  discussed  earlier.  Another  way  of  considering  this,  is  that  the  masses 
were  adjusted  until  the  invariant  value  of  forced  inclination  matched  the  inclination  of  the 
outer  disk.  The  masses  of  the  planets  in  units  of  stellar  masses  and  in  imits  of  Jupiter 
masses  are  shown  in  table  8-2 

The  models  associated  with  the  various  simulations  based  on  parameters  developed  in 
this  section  are  discussed  and  shown  in  the  following  section. 

8.3  Model  vs.  Observations 

The  most  well  studied  characteristic  of  the  Beta  Pictoris  disk  is  the  observed  warping 
of  the  disk’s  mid-plane.  It  has  been  shown  in  previous  analysis  that  the  warp  can  be 
considered  as  an  edge-on  (~  2  —  3  degrees)  view  of  a  multi-component  disk  in  which  the 
components  have  a  non-zero  mutual  inclination.  The  most  recent  observations  suggest  that 
this  multi-component  system  consists  of  three  primary  components. 

The  first  goal  of  any  model  simulating  the  Beta  Pictoris  disk  must  reproduce  the  three 
major  symmetry  planes.  Per  the  discussion  in  chapter  4,  the  major  symmetry  planes  are 
likely  to  be  associated  with  specific  angular  momentum  planes  associated  with  the  planetary 
orbits.  To  recreate  this  warping,  it  is  thus  necessary  to  invoke  a  configuration  of  planetary 
orbits  corresponding  to  the  observed  symmetry  planes  of  the  inner  disk  and  whose  total 
angular  momentum  vector  defines  the  observed  large  scale  symmetry  plane.  The  presence 
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Figure  8-14:  Simulated  disk  based  on  orbital  elements  constrained  by  the  UF  18  fim  images. 
The  values  for  the  forced  orbital  elements  are  those  shown  in  figure  8-11  and  the  values  for 
the  planetary  masses  are  those  shown  in  table  8-2. 


of  planets,  however,  must  produce  additional  effects  as  the  planets  interact  with  the  circum- 
stellar  material.  These  effects  should  produce,  at  the  very  least,  gaps  in  the  circumstellar 
disks.  The  final  model  of  the  Beta  Pictoris  disk  must  therefore  reproduce  the  major  symme¬ 
try  planes,  have  planetary  orbits  consistent  with  the  observed  gaps,  and  have  forced  orbital 
elements  that  will  shape  the  disk  in  such  a  way  as  to  reproduce  the  features  observed  in 
the  high  resolution  images.  Albeit  a  non-unique  solution,  one  planetary  configuration  that 
produces  a  model  satisfying  these  requirements  is  presented  in  this  section. 
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Figure  8-15:  The  large  scale  structure  of  the  model  disk  shown  in  figure  8-14,  but  with 
each  vertical  slice  normalized  to  the  peak  value  of  the  shce. 

those  observations.  A  plot  illustrating  the  large  scale  structure  of  the  model  disk  is  shown 
in  figure  8-15.  This  image  is  clearly  consistent  with  the  outer  warp  identified  by  Heap  et 
al.,  and,  to  the  same  extent  that  the  UF  images  are,  is  consistent  with  the  inner  plane. 

The  primary  symmetry  planes  are  significant  as  they  relate  the  inner  planets  to  the 
large  scale  structure.  The  inner  and  intermediate  symmetry  planes  reveal  the  orientation 
of  the  angular  momentum  planes  associated  with  the  planets.  The  outer  or  invariant  plane 
of  the  disk  is  associated  with  the  system’s  total  angular  momentum.  These  connections  are 
important  because  they  imply  that  a  simple  way  to  make  a  first  approximation  of  a  potential 
planetary  system  can  be  accomplished  by  simply  considering  the  angular  momentum  for 
each  planetary  orbit  in  the  proposed  system.  Therefore  a  matching  of  major  symmetry 
planes  in  the  model  is  not  a  significant  accomplishment.  To  exemplify  this  point,  a  simpler 
model  could  easily  reproduce  the  major  symmetry  planes.  An  approach  to  modeling  the 
major  symmetry  planes  based  on  these  assumptions  was  discussed  in  section  8.2.1.  In  the 
case  of  Beta  Pictoris,  planets  were  required  to  be  chosen  in  a  way  that  the  orientation  of 
their  orbits  matched  the  inner  and  intermediate  symmetry  planes.  The  masses  then  needed 
to  be  adjusted  to  provide  the  necessary  total  angular  momentum  and  hence  the  invariant 
plane  associated  with  the  total  angular  momentum  vector. 

This  implies  that  if  only  the  three  major  symmetry  planes  were  desired,  this  could  be 
accomplished  with  2  planets  -  one  for  each  of  the  inner  symmetry  planes.  The  model  can 
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Figure  8-16:  A  simple  model  of  the  distribution  of  material  in  a  system  perturbed  by  three 
planets.  The  inner  planes  are  associated  with  planes  of  the  perturbing  planets  and  the  large 
scale  invariant  plane  is  that  associated  with  the  systems  total  angular  momentum  vector. 

be  further  tuned  by  adding  additional  planets.  Such  a  toy  model  is  shown  in  figure  8—16, 
but  as  discussed  above  is  relatively  meaningless.  Fortunately,  additional  constraints  exists. 
The  presence  of  planets  will  have  additional  effects;  such  effects  are  clearly  observed  in  the 
UF  images.  These  observable  effects  include  the  gaps  and  small  scale  structure  that  are 
discussed  in  the  next  section. 

8.3.2  Brightness  Enhancements  and  Voids 

A  model  created  using  the  forced  elements  matching  the  constraints  extracted  firom 
the  observed  sections  of  disk  is  shown  in  figure  8-14.  The  model  was  created  by  simulating 
each  section  of  disk  observed  in  the  UF  images.  Each  section  of  disk  was  influenced  by  the 
distribution  of  proper  elements  discussed  in  the  previous  chapter  and  the  forced  elements 
determined  for  the  seven  planet  system  shown  in  figure  8-11.  As  a  note,  the  five  planet 
system  was  equally  capable  of  producing  the  features.  The  distribution  of  surface  area  was 
created  by  using  a  superposition  of  each  section  of  disk  where  each  section  was  weighted 
based  on  the  estimated  form  of  the  number  density  radial  distribution  given  in  equation 
7.6.  The  final  image  was  created  using  the  flux  determined  for  10  /im  particles. 

The  resulting  model  clearly  shows  the  sections  of  disk  observed  in  the  18  nm.  images 
with  the  same  orientation  relative  to  the  reference  plane.  Furthermore,  the  oflfeet  of  the 
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center  of  symmetry  with  respect  to  the  star  suggested  for  each  section  of  disk  has  also  been 
reproduced.  The  eccentricities  associated  with  each  planet  also  are  consistent  with  the  gaps 
in  the  model  image  as  well  as  the  observations. 

In  addition  to  the  gaps  created  by  accretion,  it  is  possible  for  particles  to  enter  the 
regions  near  a  planet  that  could  remove  the  particle  form  the  system.  As  was  discussed  in 
chapter  4,  because  of  the  vector  relationship  between  the  proper  and  forced  orbital  elements 
particles  will  reach  a  maximum  eccentricity  that  may  or  may  not  bring  them  close  enough 
to  a  planet  to  be  removed  and  thus  cause  depletion  for  a  particular  section  of  disk.  This 
is  especially  true  for  particles  that  are  in  regions  influenced  by  secular  perturbations.  A 
calculation  to  assess  this  effect  is  shown  in  flgure  8-17.  The  top  panel  shows  the  forced  ec¬ 
centricity  for  the  seven  planet  system  whose  orbital  elements  are  shown  in  figure  8-11.  The 
middle  panel  shows  the  intersections  of  the  regions  assumed  to  be  affected  by  the  planets 
and  the  maximum  excursion  by  the  dust  particles.  It  is  at  these  intersections  that  depletion 
is  likely  to  occur.  The  bottom  plot  shows  the  gaps  created  in  the  radial  number  density 
distribution  given  in  equation  7.6.  Further  investigation  into  the  consequences  of  this  anal¬ 
ysis  in  terms  of  disk  structure  and  stability  are  saved  for  a  later  work.  Other  contributing 
influences,  such  as  resonnance  overlap,  are  also  likely  to  play  a  significant  role  and  need  to 
be  considered. 


8.4  Discussion 

In  this  chapter,  it  has  been  shown  that  the  type  of  mid-plane  warping  observed  in  the 
Beta  Pictoris  circumstellar  disk  can  be  accounted  for  by  the  secular  perturbations  associated 
with  a  multiple  planet  system.  Mid-plane  warping  is  clearly  a  consequence  of  a  distribution 
of  forced  orbital  elements  in  which  the  forced  inclination  varies  with  semi-major  axis.  This 
is  consistent  with  both  theoretical  modeling  and  observed  effects  within  the  solar  system 
dust  bands. 

Many  additional  factors  could  be  added  to  the  model  presented  here  to  add  further 
sophistication  to  the  modeling  process.  Some  of  these  will  be  discussed  in  the  following 
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Figure  8-17:  A  plot  describing  the  potentieil  impacts  on  the  distribution  of  material  as 
a  results  of  secular  perturbations,  especially  resonances.  The  top  plot  shows  the  forced 
eccentricity  for  the  seven  planet  model.  The  middle  plot  shows  the  radii  at  which  regions 
should  be  cleared  by  planets  and  the  maximum  and  minimum  excursions  of  dust  particles. 
The  bottom  plot  shows  the  impact  of  the  cleared  regions  on  the  assumed  radial  number 
density  distribution  for  the  disk  material. 
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chapter,  however  they  do  not  discredit  the  results  and  conclusions  presented  in  this  re¬ 
search.  The  basic  conclusions  of  this  chapter  are  the  following: 

•  The  environmental  factors  inferred  from  observations  of  the  Beta  Pictoris  disk  strongly 
suggest  that  the  disk  is  dominated  by  collisional  processes.  This  implies  a  reservoir 
of  larger  bodies  that  are  subjected  to  secular  perturbations.  This  creates  a  view  in 
which  the  morphology  of  the  disk  is  dominated  by  the  dynamics  of  the  larger  bodies 
despite  the  observations  being  dominated  by  the  smaller,  shorter  lived,  dust  particles. 

•  The  age  of  the  Beta  Pictoris  disk,  based  primarily  on  the  Hipparcos  data,  provides 
strong  support  for  the  assertion  that  the  Beta  Pictoris  disk  has  had  enough  time 
to  fully  evolve  under  the  influences  of  the  secular  pertmbations  associated  with  the 
planetary  systems  used  in  the  models  presented  within  this  chapter. 

•  Secular  pertmrbations  of  a  multiple  planet  system  produce  a  distribution  of  forced 
inclinations  that  vary  with  semi-major  axis  and  thus  produce  a  warped  mid-plane. 

•  The  secondary  effects  associated  with  the  presence  of  planets  and  suggested  by  the 
observations,  are  consistent  with  the  expected  effects  of  the  planets  used  to  create  the 
large  scale  warping. 


CHAPTER  9 

FURTHER  CONSIDERATIONS 


The  model  presented  within  this  dissertation,  though  based  on  reasonable  assump¬ 
tions,  is  somewhat  idealized.  A  more  complex  model  could  be  investigated  by  considering 
a  more  sophisticated  use  of  the  disturbing  function.  Two  examples  are  briefly  discussed  in 
this  chapter  and  involve  the  consideration  of  the  effects  of  a  massive  disks  and  the  hkely 
non-spherical  shape  of  the  central  star.  Both  effects  are  easily  handled  within  the  theory 
presented  in  earlier  chapters  simply  by  adding  either  more  perturbing  masses  or  by  adding 
additional  terms  to  the  disturbing  function. 

9.1  Massive  Disks 

It  was  asserted  earlier  that  the  mass  of  the  dust  within  the  disk  of  the  Beta  Pictoris 
model  is  insignificant  compared  to  the  combined  mass  of  the  planets.  However,  this  is  based 
on  many  assumptions,  one  of  which  is  the  assumed  value  of  the  dust  to  gas  ratio  which  has 
significant  uncertainty.  Furthermore,  estimated  values  of  the  disk  mass  based  the  sub-mm 
measurements  discussed  earlier  ignore  the  contribution  from  the  more  massive  bodies  like 
asteroids. 

Therefore,  it  is  reasonable  to  consider  the  effects  on  the  forced  orbital  parameters  when 
the  mass  of  the  disk  becomes  comparable  to  that  of  the  combined  mass  of  the  planets.  Figure 
9-1  illustrates  the  result  of  one  such  calculation.  This  figure  shows  the  forced  inclination 
calculated  for  a  massless  test  particle  presumably  within  a  disk  that  is  disturbed  by  three 
massive  bodies  on  a  similar  scale  and  with  similar  values  of  inclination  as  those  used  in 
the  Beta  Pictoris  model.  The  filled  circles  of  zero-inclination  correspond  to  masses  that 
were  used  to  simulate  the  distribution  of  mass  within  the  disk  by  equally  distributing  the 
assmned  disk  mass  among  masses.  The  three  different  realizations  allowed  for  different 
values  of  total  disk  mass.  The  top  panel  shows  the  result  when  the  mass  of  the  disk  is 
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Figure  9-1:  The  forced  inclination  for  a  test  particle  perturbed  by  a  three  planet  system; 
the  planets  are  shown  by  the  filled  circles  with  non-zero  values  of  inclination.  The  mass  of 
the  disk  is  simulated  by  adding  additional  bodies  (the  zero  inclination  objects).  The  top 
panel  shows  the  result  when  the  combined  mass  of  the  disk  is  equal  to  10“®  that  of  the 
combined  mass  of  the  planets.  The  middle  panel  shows  the  same  for  a  mass  ratio  of  10“^ 
and  the  bottom  shows  the  same  for  a  mass  ratio  of  10“^. 
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arguably  insignificant  to  that  of  the  planets;  mass  ratio  of  ~  10“®.  The  middle  panel  shows 
the  same  system  but  with  the  mass  ratio  increased  to  10”^.  The  bottom  panel  shows  the 
result  when  the  mass  ratio  is  about  10“^.  As  a  reference,  this  research  assumes  an  upper 
value  of  Mdisk/Mpianets  ~  10”^  “  10“'*,  which  is  based  on  the  estimate  of  the  total  Beta 
Pictoris  disk  mass  of  ~  0.1  - 1.0  M®.  Though  this  value  neglects  the  potential  contribution 
from  asteroids  and  other  massive  bodies  within  the  disk,  a  significant  enough  margin  of 
error  exists  such  that  the  results  obtained  are  reasonably  valid. 

The  results  shown  in  figure  9—1  show  the  characteristics  that  should  be  expected  based 
on  discussions  in  previous  chapters.  First,  for  a  common  semi-major  axis  the  forced  orbital 
elements  for  the  test  particle  are  equal  to  the  osculating  orbital  elements  of  the  massive 
bodies.  Second,  the  region  around  which  the  massive  bodies  impose  their  orbital  elements 
on  the  test  particles  is  a  strong  function  of  mass.  Therefore,  as  the  mass  of  the  disk  increases, 
the  influence  of  the  disk  becomes  comparable  to  that  of  the  planets.  This  causes  the  forced 
elements  to  match  not  only  the  planets,  but  the  disk  as  well. 

The  net  results  of  the  massive  disk  is  to  dampen  the  effects  of  the  planets.  This  is 
seen  by  the  fact  that  the  regions  over  which  the  planets  exert  their  influence  becomes  more 
spatially  limited  as  the  disk  mass  increases.  Another  eflfect  is  that  the  invariant  plane  of 
the  disk,  defined  by  the  asymptotic  value  of  forced  inclination,  approaches  a  value  closer  to 
that  of  the  disk.  This  would  cause  large  scale  warping  to  become  less  obvious  in  external 
observations.  When  coupled,  these  two  effects,  dampening  of  the  structure  created  by  the 
planetary  perturbations  and  forcing  the  invariant  plane  of  the  system  and  the  plane  of  the 
disk  to  become  coplanar,  would  serve  to  diminish  the  observable  characteristics  of  planets 
existing  within  a  disk.  This  suggests  that  the  there  may  be  an  ideal  epoch  over  which 
the  observations  of  structure  within  a  circumstellair  can  be  used  to  identify  the  dynamical 
signatures  of  planets. 

Though  this  example  shows  some  of  the  important  effects  of  a  massive  disk  on  the  ob¬ 
servable  characteristics  of  a  circumstellar  disk,  it  is  still  overly  simphfied.  A  model  based  on 
additional  masses  spread  over  a  larger  region  would  better  simulate  a  realistic  distribution 
of  dust  material. 
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Figure  9-2:  An  oblate  spheroid  with  an  oblateness  of  0.25. 

9.2  Oblateness  of  the  Central  Star 

The  calculations  performed  in  this  research  are  based  on  the  assumption  that  the 
central  mass  is  perfectly  spherical.  However,  a  rotating  non-rigid  body  will  flatten  to  an 
oblate  spheroid.  This  is  seen  in  the  sun  and  the  planets  within  the  solar  system.  An  oblate 
spheroid  is  generally  defined  as  a  spherical  object  in  which  the  equatorial  radius  is  greater 
than  the  polar  radius.  An  example  is  shown  in  figure  9-2. 

The  oblateness  of  a  central  star  will  induce  changes  in  a  gravitational  potential  from 
one  based  on  a  spherical  central  object.  In  the  context  of  this  research,  the  oblateness  efiects 
are  manifested  as  an  addition  to  the  disturbing  function  (Murray  and  Dermott,  1999).  For 
particle  j,  the  part  of  the  disturbing  function  describing  the  eflfects  of  oblateness  can  be 
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expressed  as 


1  2  2 


2 


^22 

2^3 


2 


(9.1) 


where  Rc  is  the  radius  of  the  star  and  aj  is  the  semi-major  axis  of  the  particle.  The  J 
terms  are  the  gravitational  moments.  These  terms  can  be  expressed  as  functions  of  the 
oblateness,  e  =  (o  —  c)/a,  and  the  oblateness  constant,  q  =  r^u}^/{GM).  Typical  solar 
values  for  oblateness  are  thought  to  lie  between  5.0±  0.7x10“®  (Dicke  and  Goldenberg 
1967)  and  9.6±6.5  xl0“®  (Hill  and  Stebbins  1975).  Values  for  q  are  t5rpically  one  or  more 
orders  of  magnitude  smaller.  Though  Beta  Pictoris  may  have  values  significantly  different 
firom  the  sun  because  of  the  difference  in  mass  (~  I.SMq)  and  a  different  rotation  rate, 
these  values  should  give  an  idea  of  the  general  magnitude. 

Based  on  the  model  parameters  considered  in  this  research,  oblateness  should  have 
minor  effects  on  the  results.  The  J  terms  are  likely  to  be  in  the  realm  of  ~  10~®.  The 
mininniim  value  of  semi-major  axis  considered  for  particles  within  the  disk  was  ~  20  AU 
and  the  the  radius  of  Beta  Pictoris,  if  comparable  to  the  sun,  is  ~  10“®  AU.  This  implies 
that  the  the  terms  J^{Rc/a)'^  are  ~  (10“®)”(10“'*)”^.  When  compared  to  the  other  terms 
in  the  disturbing  function,  these  are  clearly  negligible. 

The  appropriate  terms  to  consider  in  making  the  comparison  are  the  elements  of  the 
A  and  B  matrices  associated  with  the  secular  part  of  the  disturbing  function. 
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where  the  values  of  A  and  B  are  functions  of  mass  and  semi-major  axis  only,  and  are  given 

by 
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(9.3a) 

II 

(9.3b) 
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These  values  can  be  thought  of  as  constant  matrix  elements  and  have  typical  values  of 
~  10~^,  for  the  type  of  configuration  of  planets  used  in  the  final  models.  These  values  are 
significantly  larger  than  those  associated  with  the  oblateness.  This  suggests  that  for  the 
conditions  and  scales  considered  in  this  research,  oblateness  plays  a  small  or  insignificant 
role. 

Though  oblateness  may  not  play  a  role  in  Beta  Pictoris,  the  effects  of  a  significantly 
oblate,  rapidly  rotating  star  on  the  accompanying  circumstellar  disk  could  prove  to  be  an 
interesting  problem  to  study. 


CHAPTER  10 
CONCLUSIONS 


The  research  presented  in  this  dissertation  was  an  investigation  into  the  phenomenon 
of  mid-plane  warping  in  circumstellar  disks.  This  particular  effort  was  inspired  by  the  ob¬ 
servations  made  by  the  UF  mid-IR  team  of  the  circumstellar  disk  of  Beta  Pictoris.  These 
observations  provided  images  that  showed  with  rmprecedented  clarity  the  featmres  of  the 
disks  mid-plane.  Within  the  images,  was  clear  evidence  of  a  complex  morphology.  Fur¬ 
thermore,  analysis  of  the  brightness  distribution  suggested  regions  in  which  the  dust  was 
depleted.  Both  of  these  characteristics  strongly  suggest  the  presence  of  planets. 

The  basic  circumstellar  environment  was  modeled  using  values  for  input  parameters 
taken  from  the  extensive  coverage  of  Beta  Pictoris  found  in  the  literature.  These  parameters 
included  the  dust  grain  material  properties,  the  radial  number  density  distribution,  and 
the  proper  orbital  elements.  Also  taken  from  the  literature  was  the  data  necessary  to 
estimate  the  timescales  associated  with  different  physical  processes  and  in  turn  validate  the 
applicabihty  of  a  secular  analysis. 

Using  secular  perturbation  theory,  a  circumstellar  disk  perturbed  by  a  system  of  planets 
was  simulated.  The  planets  considered  were  constrained  by  observable  features  of  the  disk. 
These  features  included  both  large  scale  and  small  scale  features.  The  large  scale  features 
were  the  major  symmetry  planes  of  the  system  with  transitions  at  ~20  and  80  AU  and  have 
been  documented  in  multiple  sources.  The  small  scale  features  were  clearly  visible  only  in 
the  UF  images  and  consisted  of  brightness  enhancements  and  small  scale  warping. 

The  large  scale  symmetry  planes  were  addressed  first.  Analysis  of  these  planes  provided 
insight  into  the  dominating  angular  momentmn  planes  as  well  as  the  system’s  total  angular 
momentum  plane.  The  inner  symmetry  planes  constrained  the  orientation  of  the  planetary 
orbits  and  the  outer,  invariant  plane  of  the  disk  constrained  the  relative  masses  of  the 
planets  by  constraining  the  total  system  angular  momentum. 
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The  small  scale  features  constrain  the  disk  in  two  ways.  First,  the  regions  corresponding 
to  potential  gaps  were  likely  to  be  related  to  the  location  of  the  planets.  Based  on  the 
observed  enhancements  in  brightness  and  the  interpretation  that  the  enhancements  were 
visible  sections  of  disk,  the  potential  gap  regions  were  isolated.  The  radial  width  of  these 
regions  were  an  upper  limit  to  the  expected  gap  width  associated  with  a  planet,  w  ~  re. 
These  gap  regions  also  provided  a  constraint  to  the  masses  through  the  assessment  of 
the  accretion  zones.  The  accretion  zones  are  related  to  the  Hill  sphere  associated  with 
the  corresponding  planet.  The  radius  of  the  sphere  is  determined  by  the  planetary  mass 
and  the  distance  from  the  star.  Therefore  the  gaps  also  provides  an  upper  limit  to  the 
mass  a  planet  presumed  to  be  located  in  that  particular  gap.  The  second  constraint,  was 
associated  with  the  locations  of  the  brightness  enhancements  with  respect  to  the  reference 
mid-plane.  The  inclination  and  orientation  of  a  disk  section  with  respect  to  the  mid-plane 
of  the  disk  is  related  to  the  local  value  of  the  forced  inclination  and  longitude  of  ascending 
node  respectively.  This  constraint  on  the  forced  inclination  and  node  is  a  constraint  on  the 
magnitude  of  the  perturbation  and  hence  a  constraint  on  the  inclinations  and  orientation 
of  orbits  of  the  perturbing  planets. 

This  research  suggests  that  the  quantifiable  characteristics  of  the  features  observed  in  a 
circumstellar  disk  can  be  used  to  extract  information  that  will  constrain  the  characteristics 
of  a  perturbing  planetary  system.  This  means  that  high  resolution  images  of  circumstellar 
disks  can  ofier  not  only  an  alternative  to  methods  like  stellar  radial  velocity  measurements 
as  a  means  to  detect  exosolar  planets,  but  can  actually  provide  solid  constraints  on  the 
characteristics  of  what  are  likely  to  be  very  complex  dynamical  systems. 


APPENDIX: 

PLANETARY  ORBITAL  ELEMENTS 


Throughout  this  research  comparisons  have  been  made  between  proposed  planetary 
systems  and  the  solar  system.  The  following  table  provides  the  orbital  elements  associated 
with  each  planet  in  the  solar  system  (excluding  Pluto)  (Murray  and  Dermott,  1999). 


Planet  op  (AU)  m  (lO^^kg)  eg  Ip  (deg.)  tup  (deg.)  flp  (degT) 


Mercury 

0.38709893 

0.3302 

0.20563069 

7.00487 

77.45645 

48.33167 

Venus 

0.72333199 

4.8685 

0.00677323 

3.39471 

131.53298 

76.68069 

Earth 

1.00000011 

5.9736 

0.01671022 

0.00005 

102.94719 

348.73936 

Mars 

1.52366231 

0.64185 

0.09341233 

1.85061 

336.04084 

49.57854 

Jupiter 

5.20336301 

1898.6 

0.04839266 

1.30530 

14.75385 

100.55615 

Saturn 

9.53797932 

568.46 

0.05415060 

2.48446 

92.42194 

113.71504 

Uranus 

19.19126393 

86.832 

0.04716771 

0.76986 

170.96424 

74.22988 

Neptune 

30.06896348 

102.43 

0.00858587 

1.76917 

44.97135 

131.72169 

The  above  orbital  elements  are  for  the  epoch  of  J2000  with  respect  to  the  mean  ecliptic 
and  the  equinox  of  J2000. 
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